Considerable Uncertainty of Simulated Arctic Temperature Change in the Mid‐Holocene Due To Initial Ocean Perturbation

Arctic temperature is one of the most uncertain aspects of mid‐Holocene (MH) climate change modeling, usually attributed to the different responses of different models to external forcing. However, in this study, we find that significant discrepancies (i.e., the noise is close to the signal in term of climate change) in the MH Arctic temperature changes can occur within the same model and for identical external forcing due to initial ocean condition perturbations. It is shown that initial ocean perturbations can affect the surface energy budget change through the uncertain cloud effect on shortwave radiation in boreal summer. The resulted uncertain change in summer surface heat flux alters the subsequent autumn and winter sea ice and contributes to significant differences in Arctic temperature via sea ice‐albedo feedback. This study suggests that internal uncertainty of an individual model is a non‐negligible source of overall uncertainty in simulating the MH Arctic temperature change.


Introduction
The mid-Holocene (MH) was the most recent warm geological period in Earth's history around 6,000 years ago.Proxies indicate that the global temperature was ∼0.5°C warmer than the preindustrial (PI) era, with particularly pronounced warming over the North Hemisphere (NH) high latitudes (Bartlein et al., 2011;D. Kaufman et al., 2020;Marcott et al., 2013;W. Zhang et al., 2022).Due to the spatial-sparse distribution and limited variables of proxies, the underlying dynamic processes remain elusive.Climate models have emerged as a valuable tool for paleoclimate research, which provides the quantitative signal of climate change and more importantly the physically explainable mechanism (e.g., Braconnot et al., 2012).Cross-validation between climate models and proxies not only helps to reveal the drivers for past climate change, but also tests whether the responses of climate models to various forcings are correct, determining the credibility of their projected future climate change.However, climate models generally failed to capture the warmer MH recorded by proxies (Kaufman & Broadman, 2023;Lohmann et al., 2013).For example, climate models participating in the Paleoclimate Model Intercomparison Project Phase 4 (PMIP4) performed an average of 0.3°C colder in the MH-none of them simulated a warmer Earth (Brierley et al., 2020).Such discrepancy is a part of the "Holocene temperature conundrum": proxies indicate a temperature maximum at early and middle Holocene and cooling afterward, while climate models suggested a continuous warming Holocene (Z.Liu et al., 2014).
Great efforts have been made to explain and reduce the gap between proxies and model results.The Holocene annual SST showed a warming trend by transforming the seasonal SST records to annual mean (Bova et al., 2021), showing that the puzzle is partly due to the seasonal bias in the interpretation of proxies (Z.Liu et al., 2014).Holocene temperatures were cooling in the NH high latitudes and warming in the tropics, emphasizing a latitudinal dependent response (Bader et al., 2020;G. Lohmann et al., 2013).Assimilating proxies to a climate model, the Holocene temperature was rising, and the MH was ∼0.1°C colder than PI (Osman et al., 2021).It was concluded that uncertainty is due to the seasonal bias and poor spatial coverage of the data.The model-data inconsistency could also result from unrealistic responses of external forcings or lack of forcings in climate models.After removing the MH atmospheric dust load caused by enhanced monsoon and resulting vegetation cover, the global temperature increased by 0.3°C (Y.Liu et al., 2018).Park et al. (2018Park et al. ( , 2019) ) emphasize the role of the Arctic sea ice loss in MH and its feedback on NH temperature.If climate models can simulate the mid-latitude Arctic as warm as proxies suggest, they should indicate a warmer NH than PI.Recently, Thompson et al. (2022) applied expanded vegetation in the Sahara and NH mid-latitudes during the early and middle Holocene in their model, whereupon the global temperature increased by 0.8 and 0.7°C, respectively, especially over the Arctic, which is over 5°C warmer.
Though the conundrum remains puzzling (Bova et al., 2021;F. Chen et al., 2023;Kaufman & Broadman, 2023;Osman et al., 2021), plenty of studies pointed out that the Arctic could be the key, where the polar amplification effect acts as a bridge between the external forcing and global temperature in various climate models (e.g., Bader et al., 2020;J. Chen et al., 2022;Park et al., 2018Park et al., , 2019;;Routson et al., 2019;Thompson et al., 2022).Nevertheless, responses of MH Arctic temperature changes are largely inconsistent among different models (Brierley et al., 2020), highlighting the uncertainty of models in simulating physical processes over the Arctic and emphasizing the uncertainties regarding the Arctic.
The uncertainty of a multi-model ensemble can be due to the models' different responses to external forcing, linked to their dynamical cores and physical parameterizations.Meanwhile, the uncertainty within an individual model (hereafter as internal uncertainty), for example, resulting from different initial conditions, should also contribute.Since the equilibrium simulations in PMIP were generally integrated over hundreds to thousands of years and the last hundreds of years were used for the analysis, the influence of the initial condition was considered ignorable.Nevertheless, in this study, using a set of MH and PI experiments conducted with identical model setups and perturbing only the initial ocean condition, we show that the Arctic temperature changes in the MH exhibit considerable inconsistencies.This suggests a non-negligible influence of internal uncertainty on the modeling of Arctic temperature change in the MH.

Data and Methods
The simulations are performed with the third generation of climate model developed at the Alfred Wegener Institute (AWI-CM3), consisting of the atmosphere, ocean and sea ice components.The atmosphere model is the Integrated Forecast System in the version OpenIFS (version cy43r3v2;ECMWF, 2017aECMWF, , 2017bECMWF, , 2017c) ) developed at the European Center for Medium-Range Weather Forecast.The ocean model is the second version of the Alfred Wegener Institute's global ocean model, the Finite volumE Sea ice-Ocean Model (FESOM2) with unstructured meshes, which contains the sea ice module (Danilov et al., 2017).For all the experiments, the OpenIFS has a TCo95 horizontal resolution (∼100 km) with 91 vertical layers, coupled with the FESOM2 using the CORE2 mesh (∼0.13 million surface nodes, 47 vertical layers).The AWI-CM3 has a good reproduction of observed climate in many aspects (Shi et al., 2023;Streffing et al., 2022).For more model details, we refer to https://awi-cm3-documentation.readthedocs.io/en/latest/index.html.
First, a PI spin-up experiment is integrated for 700 years under PI boundary conditions.A MH spin-up experiment, which is restarted from the 500th year of the PI spin-up ocean state is integrated for 200 years, forced by the MH boundary conditions (labeled as 501-700th simulation year; Figure S1 in Supporting Information S1).After the PI and MH spin-up experiments, we run the PI and MH experiments for 100 years under respective forcings (labeled as 701-800th simulation year), which are used for the analysis.For each period, we create 5 ensemble runs with initial perturbations based on restarts from the 690th, 695th, 698th, 699th, and 700th year ocean states, respectively.The PI and MH boundary conditions are listed in Table S1 in Supporting Information S1, with differences in orbital parameters and greenhouse gas (GHG) concentrations.The global annual mean surface air temperature (GAST) evolution shows little trends for the last 100 years, implying that the simulations reach quasiequilibrium (Figure S1 in Supporting Information S1).
All members of the AWI-CM3 ensemble are considered equally and independently, giving us 25 (5 × 5) differences between the MH and PI experiments.For convenience, these differences are named the "difference ensemble (DE)."Their arithmetic median is defined as the signal of the climate change caused by the external forcing, and their standard deviation is considered as the noise depending on the initial conditions.The (absolute) signal-to-noise-ratio (SNR) is used to measure the robustness of the modeled climate change.Since this ensemble is artificially generated, its members exhibit strong auto-correlation.To assess the significance level of our difference ensemble, we employed Monte Carlo simulation.Initially, we generated two MH-PI difference ensembles, each consisting of 25 values obtained by subtracting 5 random MH values from 5 random PI values.We then calculated the correlation coefficient between these two 25 values and repeated these steps 10,000 times.Subsequently, we determined the 95th percentile of the 10,000 correlation coefficients, which served as the threshold for conducting a 95% significance test when performing correlation analyses within the ensemble.

Uncertainty in Modeling the Mid-Holocene Temperature
The GAST of the AWI-CM3 ensemble ranges from 12.86 to 12.93°C for MH and 13.15-13.22°Cfor PI, implying an average cooling of 0.29°C in MH relative to PI (Figure S2 in Supporting Information S1).This value is close to the PMIP4 ensemble mean ( 0.3°C; Tian et al., 2022).The spatial pattern of surface temperature change exhibits a large-scale cooling over the tropical and NH mid-latitudes and a warming over NH high latitudes (Figure 1a), resembling the PMIP2, PMIP3, PMIP4 ensemble result (Brierley et al., 2020;G. Lohmann et al., 2013).The general global cooling is mainly caused by the lower GHG concentrations in MH.The most significant cooling takes place over the North African and South Asian monsoon regions, which is mainly caused by the increased cloud induced by the enhanced monsoonal precipitation in response to stronger summer solar insolation (Figure S3 in Supporting Information S1).The mechanisms can be attributed to the insolation forcing in interglacials (Herold & Lohmann, 2009).The warming over the Arctic is caused by sea ice loss due to increased summer insolation and its feedback to the atmosphere.Our seasonal temperature changes generally agree with the PMIP4 results (Figure S4 in Supporting Information S1), following the seasonal solar insolation, suggesting that AWI-CM3 largely reproduces the MH temperature responses similar to other models.
However, the GAST changes between the MH and PI differ notably for the AWI-CM3 initial ocean perturbation experiments, ranging from 0.37 to 0.22°C among the DE (Figure S2 in Supporting Information S1).The GAST difference between the warmest and coldest member (0.15°C) exceeds one standard deviation among the PMIP4 ensemble (0.12°C; Brierley et al., 2020;Tian et al., 2022).This suggests that MH climate change uncertainty also exists within a single model, which could lead to large uncertainty in a multi-model ensemble result because most models only provide one ensemble for the multi-model comparison.
Specifically, the noise of GAST changes in the DE is 0.04°C and the corresponding SNR is 7.3, which implies that the global mean temperature response among AWI-CM3 simulations is less affected by the initial ocean perturbation.Nevertheless, this large global SNR is mainly dominated by the consistent strong cooling over the tropics, especially over the NH monsoon regions and tropical eastern Pacific, where the local SNR exceeds 20, indicating the robustness of AWI-CM3 in simulating the MH climate change over these regions (Figure 1).Meanwhile, the SNR reduces sharply from tropical to polar regions in both hemispheres, with SNR close to or even less than one over the Arctic/Antarctic, emphasizing the pronounced influence of the initial ocean on the simulated MH polar climate change.For the seasonal temperature change, the SNR shows similar spatial distributions (Figure S4 in Supporting Information S1).The smallest SNR over the Arctic (Antarctic) occurs in local winter, which implies that this large uncertainty of polar temperature response may be closely related to the seasonal sea ice variation.Although both the Arctic and Antarctic show widespread temperature changes in the DE, the global temperature change is dominated by the Arctic (Figures S1 and S2 in Supporting Information S1).The inter-DE Arctic temperature change is approximately 6 times of the global mean temperature change (Figure S5 in Supporting Information S1).This ratio is significantly larger than the well-known Arctic amplification effect, which is about twice as warm as the global average (Winton, 2006), indicating additional processes happening there other than the cooling effect due to reduced GHG concentrations (Brierley et al., 2020).In Antarctica, on the other hand, the warming/cooling rate is not significantly different from the global average.Its zonal mean is close to 0, contributing little to the global temperature difference in the DE.

Arctic Sea Ice Responses
To better illustrate the different temperature responses due to the initial perturbation, the five coldest and five warmest ensemble members, regarding GAST, of DE are chosen (Figure 2).Generally, the seasonal cycle of MH temperature changes directly follows the solar insolation because the GHG effect is more spatial-and seasonal-

Geophysical Research Letters
10.1029/2023GL106337 uniform (Figure S3 in Supporting Information S1).For the Arctic region, the largest warming occurs in boreal autumn, lagging the orbital-induced solar insolation change by 3 months.This is because the increased solar insolation will continuously cause an Arctic sea ice loss during boreal summer, hence less sea ice accumulation in autumn, then strongly warming the Arctic atmosphere (e.g., Park et al., 2019).For the coldest ensemble, the autumn warming lasts until December, and the winter is cooling again by the reduced solar insolation and accompanied sea ice increase (Figure 2b).However, the warmest ensemble exhibits a long-lasting warming from autumn to February, and the winter cooling is also less remarkable (Figure 2a).As a result, the warmest ensemble displays significantly warmer NH high latitudes throughout the year (Figure 2c).
Because changes in Arctic sea ice play an important role in altering the temperature response to solar insolation, we also show the Arctic sea ice concentration (SIC) in the DE (Figure 2d).The Arctic sea ice loss in MH shows a large spread, especially during boreal autumn.In the coldest ensemble, the Arctic SIC reduced by 7% in September, which is one-third less compared to a 10.5% SIC reduction in the warmest ensemble.This is directly reflected in their autumn temperature responses (Figure 2c).Moreover, in the coldest ensemble, the Arctic SIC stops reducing around December and turns to increasing from January to May.This sea ice growth in the latter period is responsible for the colder Arctic in winter and spring.The SIC in the warmest ensemble decreases throughout the year, which is distinct form that in the coldest ensemble.Particularly, the Arctic sea ice remains reduced in winter despite the solar insolation decreases, and decreases significantly in the sequent summer, warming the Arctic.These differences in seasonal SIC changes are also obvious in the spatial distribution (Figure 3).Compared to the coldest ensemble, the warmest ensemble simulates a more pronounced sea ice loss in September, concentrating on the sea ice edge, where the sea ice is thinner and vulnerable to solar insolation change (Figures 3a-3c).In March, the climatological sea ice extends southward.The coldest ensemble indicates sea ice loss mainly over the Barents Sea, which is insignificant in the warmest ensemble.The sea ice also loses remarkably over the Bering Sea in both coldest and warmest ensembles in response to the decreased solar radiation through winter (Figures 3d and 3e).However, along with the seasonal movement of sea ice, the increased loss of summer sea ice in the warmest ensemble is unfavorable for winter sea ice growth over the Bering Sea (Figure 3f).In summary, different sea ice loss responses to summer solar insolation change between MH and PI accounts for the uncertainty in simulating Arctic temperatures change.

Surface Energy Budget
An unresolved question remains: why do MH changes in sea ice and temperature in the Arctic differ significantly between simulations using the same model and identical external forcing with different initial conditions?We analyze the surface energy budget over the Arctic, which is the direct controller of sea ice amount.In the individual PI and MH simulations, the climatological surface heat flux is quite consistent for each component (Figures 4a and 4b): the net shortwave radiation (NSR), net longwave radiation (NLR), sensible heat and latent heat at the surface (Equation 1).
Heat Flux surface = NSR surface + NLR surface + Sensible Heat + Latent Heat (1) Among these four items, the surface NSR dominates the annual cycle of heat flux.In the MH, the surface heat flux is more negative in summer and more positive in winter, indicating more surface energy gained and lost respectively relative to the PI (Figure 4c).Besides that, the DE exhibits considerable noise compared to the signal in the heat flux, particularly during summer and winter.The first noise is inherited from the NSR change while the second noise is mainly caused by the NLR and Sensible Heat changes.As analyzed above, since summer solar insolation is crucial to the simultaneous and delayed sea ice changes, we focus on the uncertainty in the summer surface NSR change.The surface NSR consists of two items: the downwelling shortwave radiation (DSR) and the upwelling shortwave radiation (USR) at the surface (Equation 2).The surface DSR and USR changes are anti-correlated since the USR is mainly the reflection of DSR by Arctic sea ice (Figure 4d).In late summer, this anti-correlation is disturbed because the NLR is less positive, implying that the atmosphere warms the surface and causes additional sea ice loss (Figure 4c).Generally, the surface NSR change is mostly triggered by the surface DSR change.However, the uncertainty of surface DSR change is smaller than the NSR change, which demonstrates that the sea ice-albedo positive feedback amplifies the uncertainty from DSR.
To connect the surface DSR with external forcing, we further decompose the surface DSR into three items: (a) the DSR at the top-of-atmosphere (TOA), which is the ultimate energy source of the whole climate system; (b) and (c): the effects of atmosphere (without clouds) and cloud on DSR before it reaches the Arctic surface, respectively (Equation 3; See Text S1 in Supporting Information S1 for the definition of E atmo.and E cloud ).First, the changes in DSR at TOA are almost identical in the DE (Figure 4e).This is reasonable since this item is merely forced by the Earth's orbital changes, which are fixed in our experiments.The atmosphere effect items also show large consistency in the entire year.As a result, the uncertainty of surface DSR change is almost completely caused by noise from changes in cloud effect, particularly in summer.In short, the cloud effect is sensitive to the initial condition, which makes the solar insolation reaching the surface different in the ensemble members, although the changes in solar insolation at TOA are identical.These differences in incoming surface solar insolation are further amplified

Geophysical Research Letters
10.1029/2023GL106337 by the sea ice-albedo feedback, leading to considerable uncertainty in the Arctic surface temperature changes in the MH.
To sum up, it is shown that even for the same climate model using identical external forcing for the MH and PI, the temperature responses in the Arctic can vary significantly due to different initial ocean conditions, which is attributed to the indeterminate cloud effect on the incoming solar insolation.

Discussion and Conclusions
The cloud effect is one of the most complex and uncertain parts of the climate system (Gettelman & Sherwood, 2016;Mülmenstädt et al., 2021;Vial et al., 2013;Zelinka et al., 2020).In addition to the sea ice-albedo feedback, the cloud feedback is an important process over the Arctic, whose sign is determined by the balance of shortwave cooling and longwave heating.Although we do not analyze the detailed cloud feedback, the differential cloud effect on shortwave radiation in DE will further alter the sign or strength of the Arctic cloud feedback, which in turn will affect the Arctic climate.For instance, cloud liquid particles are smaller than cloud ice particles and more effective at diffusing and reflecting solar insolation (Goosse et al., 2018).Due to the increased solar insolation, the warmer MH summer may contain more mixed-phase clouds and the proportion of cloud liquid particles increases as well, which will increase the planetary albedo and cool the Arctic (Bodas-Salcedo et al., 2016;Mitchell et al., 1989).These cloud feedbacks make their effects on Arctic temperature quite complicated and increase the uncertainty in the simulation.
The large uncertainty in Arctic temperature change in MH is a prominent feature in multi-model comparisons: For the PMIP4 MH experiments, the standard deviation between the models is about twice the median of the models regarding the Arctic surface temperature (Figure S6 and Table S2 in Supporting Information S1).This overall uncertainty could be due to the different model structures and internal uncertainty (e.g., dependence on initial conditions), although the latter receives less attention.The internal uncertainty is difficult to measure because most models did not have ensembles for the MH and PI.Here, we provide a rough estimate of its relative contribution.As shown in Figure S6 in Supporting Information S1, the internal uncertainty could account for up to 40% of the overall uncertainty for local annual temperature in the Arctic, which could be up to over 100% for boreal winter.This indicates that the internal noise can have a significant impact on multi-model results, highlighting the importance of detecting or reducing internal uncertainty.
In conclusion, we have shown that initial ocean condition could lead to considerable uncertainty in modeling the Holocene Arctic climate change.Previous studies also revealed that the cloud effect is sensitive to the initial response to the ocean perturbation (Andrews et al., 2012;Bony et al., 2006;Vial et al., 2013;Zelinka et al., 2012), or model resolution (Lohmann et al., 2021;Shi et al., 2020), suggesting that the uncertainty in Arctic temperature responses is a common feature modeling climate change scenarios.As all the perturbation experiments run for 100 years, it is plausible that multi-centennial variability, as observed in reconstructions, may also contribute to the uncertainty in 100-year climatologies.In our simulations, we have not identified significant multi-centennial variability in Arctic temperature, suggesting that multi-centennial variability may not play a decisive role in our specific problem (Figure S7 in Supporting Information S1).However, it's worth noting that in some climate models (e.g., EC-Earth3-LR; Q. Zhang et al., 2021), multi-centennial variability is dominant in the climate system.The inability of our climate model to capture the multi-centennial variability recorded in reconstructions makes it challenging to investigate its origin and sensitivity to initial conditions.Furthermore, it is questionable whether a longer simulation is possible to reduce the influence of initial perturbation.We further show the time evolution of related shortwave radiation processes (Figure S8 in Supporting Information S1).The DSR change at TOA is almost constant for the entire simulation period, and the change in atmosphere effect on DSR also reaches quasi-equilibrium after approximately 20 years.In contrast, the change in cloud effects declines in the first 20 years and recovers to a relatively stable state after 50 years of the integration.However, its uncertainty does not show a reducing trend along with the integration length, which indicates that a longer simulation would not expect to reduce uncertainty in the MH Arctic temperature change.Since the ensemble result becomes stable after ∼50 years, we suggest that conducting an ensemble of relatively short runs with different initial ocean conditions instead of a long run for MH and PI could help qualitatively distinguish the internal uncertainty and provide a better understanding of the results between different models.

•
The Arctic temperature changes in the mid-Holocene (MH) exhibit considerable discrepancy among the experiments with initial ocean perturbation • The summer cloud effect on solar insolation plays a crucial role in modifying the Arctic sea ice change and hence the temperature response • The internal uncertainty within a single model is non-negligible in simulating the MH Arctic temperature change Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.(a) Signal, (b) noise, and (c) signal-to-noise-ratio of annual mean surface air temperature changes in the DE.Units for (a) and (b) are °C.Right panel of each subplot is the corresponding zonal mean.

Figure 2 .
Figure 2. Zonal mean surface temperature change in the (a) warmest ensemble and (b) coldest ensemble and (c) their difference (in °C) in the DE.(d) Regionally averaged (70°-90°N) sea ice concentration (SIC) changes (in %) in the warmest ensemble (red curve) and coldest ensemble (blue curve).The black contours in (a) and (b) represent zero contours.The black dashed curve in (d) represents the median SIC change and the gray shading represents one standard deviation uncertainty of SIC change in the DE.

Figure 3 .
Figure 3. September sea ice concentration (SIC) changes (in %) in (a) warmest ensemble, (b) coldest ensemble and (c) their difference in the DE.Panels (d-f) are the same as (a-c), respectively, but for March.The black contours represent the climatological sea ice edge (15% SIC) for the corresponding months of the PI ensemble median.

Figure 4 .
Figure 4. Arctic (70°-90°N) surface energy budget for (a) PI and (b) MH, and (c) their differences (MH minus PI); (d) Differences in surface downwelling (downwelling shortwave radiation (DSR)) and upwelling shortwave radiation; (e) Differences in DSR at the top of atmosphere, effect of atmosphere (E atmo. ) and cloud (E cloud ).Shadings represent one standard deviation uncertainty.Positive (negative) values in (a) and (b) represent upwelling (downwelling) radiation, sensible heat or latent heat.
Geophysical Research Letters 10.1029/2023GL106337This work was supported by the National Natural Science Foundation of China (Grants 42105046 and 42088101).The work is furthermore supported through the Helmholtz program "Changing Earth-Sustaining our Future" and the project "Abrupt Climate Shifts and Extremes over Eurasia In Response to Arctic Sea Ice Change" funded through the Bundesministerium für Bildung und Forschung (BMBF, Federal Ministry of Education and Research, Germany).Jian Shi is funded by the China Scholarship Council.Open Access funding enabled and organized by Projekt DEAL.