Mixed layer depth in the PMIP4 midHolocene simulations: comparison to proxy data in North Atlantic deep convection regions

The ocean mixed layer plays an essential role in the climate system, regulating energy fluxes at the ocean-atmosphere interface. Its representation in climate models is thus critical. Here, we evaluate the mixed layer depth (MLD) in 15 models from the Paleoclimate Modelling Intercomparison Project 4 (PMIP4) against dinocyst-based MLD reconstructions from the subpolar North Atlantic for the mid-Holocene (MH, 6 ka BP). We find a large spread in MLD responses to MH forcings across the models in the present-day deep-water formation areas, underscoring the importance of model uncertainty. Most models fail to capture the direction of MLD change, and the ensemble mean does not necessarily outperform individual models. While the ensemble mean aligns broadly with proxy data in the Nordic Seas, pronounced proxy-model discrepancy in the Labrador Sea suggests that meltwater forcing is a missing parameter, and that deep-water formation in the Labrador Sea may be particularly vulnerable under a future scenario of global warming and ice-sheet melting.

1 Introduction

The ocean mixed layer is a highly homogenised surface layer of the ocean. It is the interface between the ocean interior and the atmosphere (or sea ice), where exchanges of heat, momentum and mass take place. The mixed layer is therefore a key component of the climate system. In the subpolar North Atlantic, deep-water formation in winter is associated with a deep mixed layer created by vigorous convection that feeds the lower limb of the Atlantic Meridional Overturning Circulation (AMOC) (Lazier et al., 2001; Marshall and Schott, 1999). Hence, the variability of the mixed layer depth (MLD) in the subpolar North Atlantic has attracted much interest. Evaluating past changes of MLD provides valuable insights into the sensitivity and resilience of the climate system, as the mixed layer controls the ocean's uptake of heat and carbon and hence plays a significant role in regulating global climate (Bopp et al., 2015; Kostov et al., 2014). It also informs on how the upper ocean stratification and conditions for deep-water formation may change (Roch et al., 2023; Schulze et al., 2016).

The Paleoclimate Modelling Intercomparison Project (PMIP) aims at understanding the response of the climate system to past climate forcings that differ from the present and provides an out-of-sample test for state-of-the-art numerical models by comparing them against paleoclimate reconstructions (Kageyama et al., 2018). The PMIP simulations are a means to examine the ability of models to reproduce changes in MLD under different climate states. The Middle Holocene interval (8.2–4.2 ka BP as formally defined in Walker et al., 2019), for example, is a period of the Earth's recent past, well documented by proxy records. It is one of the two interglacial intervals included in the PMIP simulations and has been used to assess the reliability of future projections (e.g., Hargreaves et al., 2013; Harrison et al., 2014, 2015). As part of the sixth phase of the Coupled Model Intercomparison Project (CMIP6), the PMIP4 midHolocene experiment is designed to examine the impact of orbital changes at 6 ka BP, when continental configurations were almost identical to modern ones and greenhouse gas levels were comparable to preindustrial ones (Otto-Bliesner et al., 2017). Compared to the previous PMIP3-CMIP5 midHolocene experiments, which only imposed the change in orbital configuration, the PMIP4 protocol includes a realistic specification of atmospheric greenhouse gas concentrations (Otto-Bliesner et al., 2017). Additionally, the PMIP4 simulations represent a new generation of climate models with higher complexity (e.g., dynamic vegetation, ocean biogeochemistry), improved parameterisations and resolution (Brierley et al., 2020; Xu et al., 2025). The Middle Holocene orbital changes lead to a different latitudinal and seasonal distribution of solar radiation, with higher (lower) insolation in the Northern Hemisphere summer (winter). The lower greenhouse gas concentrations would lead to slightly colder climates globally. In response to these external forcings, models simulate some global climate changes that are broadly consistent with theory and reconstructions, such as warmer Northern Hemisphere summer and associated shifts in tropical rainfall, reduced Arctic summer sea ice extent, cooler annual mean temperature globally but warmer in the Arctic, etc. (Brierley et al., 2020). Yet, the changes in MLD remain unexplored.

While PMIP equilibrium simulations provide a benchmark for mean-state changes, they may not accurately capture important feedbacks due to uncertainties in the boundary and initial conditions (Otto-Bliesner et al., 2017). The known uncertainties include, for instance, the prescription of modern vegetation cover and ice sheets, in addition to the lack of transient history of freshwater release from ice sheets into the ocean. Marine proxy records have demonstrated a significant impact of ice-sheet meltwater on ocean circulation during the last deglaciation (e.g., Alley and Clark, 1999; and references therein). Various model-based studies of the modern ocean also show that increasing meltwater from Greenland can lead to significant slowing of the AMOC (e.g., Golledge et al., 2019) and has the potential to disrupt deep-water formation over extended time scales (Hoshyar et al., 2025). In particular, Pontes and Menviel (2024) show that improved data-model agreement on AMOC variations over the historical period is achieved when an upper-end meltwater input is included in the simulations.

Past changes in deep convection have been difficult to reconstruct as the existing proxies are not direct indicators of convection. Recently, Wu et al. (2025) suggested that dinoflagellate cyst (dinocyst) assemblages can be influenced by wintertime deepening of the MLD in response to deep convection events, allowing for the quantitative reconstruction of MLD from dinocyst records recovered from marine sediment cores. The results from the study by Wu et al. (2025) indicate a westward migration of deep convection centres in the subpolar North Atlantic around 6 ka BP, and a late onset of modern-like deep convection in the Labrador Sea about 4 ka BP. This is compatible with previous studies using bottom water proxies such as sortable silt and benthic foraminifera to reconstruct the overflow of deep water from the Nordic Seas via the Iceland-Scotland Ridge, which indicate that deep-water formation rate in the Nordic Seas reached a maximum at 6–7 ka BP before decreasing progressively (e.g., Kissel et al., 2013; Rasmussen et al., 2002; Thornalley et al., 2013). Other studies, such as those of Renssen et al. (2005) and Thornalley et al. (2013), used transient simulations with individual models forced by Holocene changes in orbital parameters, Greenhouse gas concentrations, and ice-sheet meltwater to simulate changes in winter convection depth in the Nordic Seas and the Labrador Sea. These studies also suggest a decrease in deep convection in the Nordic Seas and an opposite trend in the Labrador Sea during the Holocene.

Here, for the first time, we compare quantitative reconstructions of winter MLD in the subpolar North Atlantic with PMIP4 model simulations for the mid-Holocene (MH, defined here as the time slice centred on 6 ka BP). We examine the differences between 15 PMIP4 models and use a combination of qualitative and quantitative approaches to assess the proxy-model agreement. We also include 13 models from the previous PMIP3 ensemble to allow for a comprehensive assessment of model robustness across different model generations. We supplement our analysis with transient simulations of the last deglaciation (PMIP4 LDv1; Ivanovic et al., 2016) to assess the impact of ice-sheet meltwater feedback.

2 Materials and Methods

We used outputs from the pre-industrial control experiment (piControl; PI) and the mid-Holocene experiment (midHolocene; MH) of 15 global climate models participating in the PMIP4 experiments (Table 1) and 13 models for the equivalent PMIP3 experiments (Table 2). The protocol and experimental design are described by Eyring et al. (2016) and Otto-Bliesner et al. (2017) for the PMIP4, and by Taylor et al. (2012) and Braconnot et al. (2011) for the PMIP3. The PI experiment represents 1850 CE and is the reference simulation from which the MH experiment branches and against which it is compared. The MH experiment is run with orbital parameters and greenhouse gas concentrations for 6 ka BP; the rest of the forcings and boundary conditions, namely the solar constant, paleogeography, ice sheets, vegetation, and aerosols, are the same as PI. The models' online calculations of MLD were not used in this study because their definitions of MLD differ. We recomputed the MLD with a potential density threshold of 0.03 kg m⁻³ and a reference depth of 10 m, as proposed by de Boyer Montégut et al. (2004) and recommended by Treguier et al. (2023) for model intercomparisons as a cost-efficient and relatively accurate estimation of the MLD. Monthly outputs of temperature and salinity from the last 100 years of each simulation were used, unless fewer years were available (PMIP3 HadGEM2-CC, 35 years available for the MH). No calendar adjustments were applied because they could affect the accuracy of variables that change abruptly throughout the year rather than gradually (Brierley et al., 2020).

Table 1 Summary of PMIP4 model configurations and features.

Table 2 Summary of PMIP3 model configuration and features.

We also included the PMIP4 last deglaciation experiment (LDv1) for process-based validation related to transient meltwater forcing. The experimental design and boundary conditions are described in Ivanovic et al. (2016). A full analysis of the participating models is presented in Snoll et al. (2024a). We only analysed simulations covering the MH in this study (Table 3). The simulations are prescribed with time-varying orbital forcing, greenhouse gases, ice sheets and other geographical changes. The GLAC-1D (Ivanovic et al., 2016) and ICE-6G_C (Peltier et al., 2015) are the two recommended ice-sheet reconstructions in the PMIP4 last deglaciation protocol, although the older ICE-5G (Peltier, 2004) is also used by some models. The protocol recommends calculating the meltwater flux based on ice volume change from the chosen ice-sheet reconstruction and applying this flux either as a globally uniform freshwater flux (melt-uniform) or at individual meltwater river outlets along the coast (melt-routed). However, the meltwater input used in the TraCE-21ka simulation does not follow the ice-sheet reconstructions mentioned above; in this simulation, a freshwater flux is used to regulate the AMOC during the early last deglaciation, so that the model produces climate changes as recorded by proxies (Liu et al., 2009; Snoll et al., 2024a). The three versions of the MPI-ESM-CR model, P1, P2 and P3, correspond to different model physics including tuning of parameters, parameterisations, and bug fixes (Kapsch et al., 2022). The output available from these long transient simulations is sometimes limited to annual and decadal averages. Recomputing the MLD from annual or decadal temperature and salinity data would probably introduce more bias than benefits since MLD is a nonlinear function of the density profile (e.g., de Boyer Montégut et al., 2004; Treguier et al., 2023). Therefore, we used the LDv1 models' online-calculated MLDs, even though their definitions differ.

Table 3 Summary of PMIP4 last deglaciation model configurations.

The MLD and sea ice reconstructions are based on 24 Holocene records of dinoflagellate cyst (dinocyst) assemblage data, compiled and quality-checked by Wu et al. (2025), from modern deep convection regions in the subpolar North Atlantic. The reconstructions were run using the modern analogue technique and the analogue R package (v0.18.0; Simpson and Oksanen, 2025). The modern analogue technique is based on the similarity between fossil assemblages and modern assemblages, using the closest modern analogues to a fossil sample to predict the corresponding environmental parameters of the past. For consistency, we used a calibration dataset extracted from de Boyer Montégut (2023), which is a monthly climatology of MLD based on temperature and salinity profile data from years 1970 to 2021, with MLD values computed using the same criterion as the one we used to compute MLD from model outputs (i.e. de Boyer Montégut et al., 2004). The modern calibration data of dinocyst species composition and sea ice concentration are from de Vernal et al. (2020). Based on cross-validation tests using the leave-one-out method, the root mean square error of prediction (RMSEP) for the reconstructions is estimated to be 40.9 m for winter MLD, 16.0 m for annual MLD, and 17.6 % for winter sea ice concentration.

To reduce the effects of systematic biases in numerical models and proxy reconstructions, results are presented as MH minus PI anomalies for proxy-model comparisons. For the proxy records, we use the average of the interval 6.5–5.5 ka BP to represent the MH; the precise signal of the PI is hard to capture given the limitations of either temporal resolution or core-top recovery, hence we use the values reconstructed from nearby surface sediment samples, which generally represent the average over the last decades to centuries. For the transient deglaciation simulations, we use the same 6.5–5.5 ka BP window to represent the MH as for the proxy records, and the average of 0.1–0 ka BP as the PI baseline value.

To quantitatively assess the proxy-model agreement, we use two complementary statistical metrics: root mean square error (RMSE) and Cohen's κ with binary categories (“positive anomaly” vs. “negative anomaly”). The proxy sites are projected to the nearest grid cell in each model's native grid prior to the analyses. RMSE is expressed in the units of the environmental variable and is therefore directly interpretable and comparable across models; however, it is dominated by amplitude errors, meaning that a model simulating little change can score better than one that reproduces the correct spatial pattern but excessive amplitude. Cohen's κ with binary categories instead evaluates only the direction of change rather than amplitude, while also corrects for the agreement expected by chance. Its limitations are that sites falling close to the category threshold are weighted the same as unambiguous ones, and that its robustness is dependent on sample size. Therefore, combining these two metrics is usually more informative than either metric alone. Additional figures and tables supporting the analyses presented here are provided in the Supplement.

3 Results

3.1 Simulated winter MLD response to midHolocene forcings

Given our primary interest, which is the implications of MLD in the deep-water formation in the MH North Atlantic, we focus here on changes in MLD during winter (January–March). The MH simulations from the PMIP4 and PMIP3 models show markedly different changes in winter MLD relative to the PI, both in terms of spatial pattern and amplitude (Figs. 1 and 2). Among the PMIP4 models, ACCESS-ESM1.5 shows the largest winter MLD anomalies, whereas INM-CM4-8 displays the smallest anomalies, especially in the North Atlantic. Compared to the PMIP4 ensemble, more models from the PMIP3 ensemble show widespread anomalies of very large magnitude (over 500 m), such as FGOALS-s2, GISS-E2-R, and KCM1.2.2.

Figure 1 Winter MLD anomalies (MH-PI) as the average of January to March in the North Atlantic from the PMIP4 ensemble (background) versus proxy data (dots).

Figure 2 Winter MLD anomalies (MH-PI) as the average of January to March in the North Atlantic from the PMIP3 ensemble (background) versus proxy data (dots).

In the North Atlantic, a recurring pattern is revealed in ACCESS-ESM1.5, EC-Earth3-Veg, GISS-E2.1G, MPI-ESM1.2-LR, and MRI-ESM2.0 from the PMIP4 ensemble: increased winter MLD in the Labrador Sea (and sometimes the Irminger Sea) concurrent with a northwestward shift in the location of MLD deepening in the Nordic Seas. The increase in MLD in the Labrador and Irminger Seas is also seen in several other models, such as CESM2, FGOALS-g3, and NorESM2-LM; on the contrary, IPSL-CM6A-LR and NESM3 are the only models that show a pronounced decrease in MLD in this area. The northwestward shift in the MLD deepening area in the Nordic Seas is also found in EC-Earth3, EC-Earth3-LR, MIROC-ES2L, and NESM3; however, FGOALS-f3-L and IPSL-CM6A-LR instead show a shift in the opposite direction. Finally, CESM2 and NorESM2-LM seem to exhibit a dipole pattern, with increased MLD in the Labrador and Irminger Seas but decreased MLD in some areas of the Nordic Seas. The anomaly pattern of increased MLD in the Labrador Sea and northwestward shift in the Nordic Seas is also found in the PMIP3 ensemble (CNRM-CM5, HadGEM2-CC, HadGEM2-ES, MPI-ESM-P).

We further analysed the seasonal cycle of MLD averaged over three deep-water formation regions as defined by Treguier et al. (2023): the Labrador Sea, Irminger Sea, and Greenland Sea (see also Fig. S1 in the Supplement). Most models simulate a seasonal cycle of the MLD comparable to that observed today: the MLD slowly deepens starting from late autumn and rapidly shoals in spring (Figs. 3 and 4). However, there is a large spread in simulated MLD values across the models in both MH and PI, with some models producing excessively deep mixed layers in some regions. In the three regions where deep convection occurs, the models seem to have their own “preference”: notable wintertime MLD deepening would generally be simulated in one or two of the three regions. In the PMIP4 ensemble, NorESM2-LM is the only model to show significant winter MLD deepening across all three regions. The other exception is EC-Earth3-LR, from which MLD deepening is not simulated in any of the three regions in MH or PI, but rather in the area extending across the Norwegian Sea, the eastern North Atlantic shelf, and the Iceland Basin (see Fig. 1). Compared to the PMIP4 ensemble, more models in the PMIP3 ensemble show significant winter MLD deepening in all three regions and all models simulate a deeper winter MLD in at least one of the regions.

Figure 3 Seasonal cycle of MLD in three deep-water formation regions from the modern climatology of de Boyer Montégut (2023) and PMIP4 model simulations of MH (a–c) and PI (d–f). The reconstructed maximums in the MH panel correspond to the maximum winter (January–March) MLD reconstructed over the MH interval in nearby dinocyst records, where available. For a map of the three regions as defined by Treguier et al. (2023), see Fig. S1.

Figure 4 Seasonal cycle of MLD in three deep-water formation regions from the modern climatology of de Boyer Montégut (2023) and PMIP3 model simulations of MH (a–c) and PI (d–f). The reconstructed maximums in the MH panel correspond to the maximum winter (January–March) MLD reconstructed over the MH interval in nearby dinocyst records, where available. For a map of the three regions as defined by Treguier et al. (2023), see Fig. S1.

3.2 Reconstructed MLD and comparison with MH simulations

The winter MLD reconstruction from dinocyst records, despite their sporadic distribution, reveals some spatially coherent patterns (Fig. 5a, c). All four sites in the Labrador Sea show shallower winter MLD during the MH. In the Nordic Seas, increased winter MLD is reconstructed at most sites from the west and southeast, while decreased winter MLD is reconstructed at the sites from the east. The PMIP4 ensemble mean MLD anomaly in the North Atlantic largely captures the recurring pattern of increase in the Labrador Sea and northwestward shift in the Nordic Seas found in individual models, and aligns broadly with reconstructions in the Nordic Seas (Fig. 5a). However, in the Labrador Sea, both the PMIP3 and PMIP4 ensembles show contrasting signals of MH anomalies compared to the reconstructions.

Figure 5 Winter (January–March) MLD anomalies (MH-PI) from the PMIP4 (a) and PMIP3 (c) model ensemble mean (background) versus proxy data (dots). Note that two of the four sites in the Labrador Sea are located very close to each other and hence superimposed in the figure. The inter-model spread is defined as the standard deviation across the PMIP4 (b) and PMIP3 (d) ensembles. (e, f) Winter (January–March) sea ice anomalies (MH-PI) from the PMIP4 ensemble mean versus proxy data and the inter-model spread.

As models often simulate MLD anomalies of the opposite sign to the reconstructions, differences in RMSE between models largely reflect differences in simulated amplitude rather than in agreement with the reconstructions. Consequently, a low RMSE may reflect a weak simulated response rather than better agreement (Sect. 2). We therefore do not use RMSE to rank model performance. The RMSE remains informative in absolute terms: for every model, it exceeds the RMSEP of the reconstruction itself, indicating that the proxy-model disagreement is larger than can be attributed to reconstruction uncertainty alone. The standard deviation across the ensemble is comparable to the ensemble-mean signal of the MLD anomaly for PMIP4, but larger for PMIP3 (Fig. 5b, d).

The PMIP4 model simulations' RMSE relative to proxy reconstructions ranges from 75.79 m of EC-Earth3 to 355.66 m of FGOALS-g3 (see Fig. 1). The ensemble mean field has an RMSE of 100.58 m. Testing the models' ability to predict the direction of change, MRI-ESM2.0 is the only model that shows statistically significant agreement with the proxy data, with a Cohen's κ coefficient of 0.44 and a 95 % confidence interval of 0.15–0.73. For the other models and the ensemble mean, the κ coefficient is generally low or negative, suggesting slight or no agreement between the proxy and the models. The confidence interval for these models and the ensemble mean includes zero (see Table S1), implying that the true κ coefficient could be zero and that the true agreement could be the same as would be expected by chance despite a positive κ value; hence, the calculated degree of agreement is not statistically significant. Among the models, ACCESS-ESM1.5 (κ=0.35) and GISS-E2.1G (κ=0.27) have the highest potential proxy-model agreement; however, this agreement may not be significant as discussed above.

None of the PMIP3 models or the ensemble mean has a statistically significant agreement with the proxy data, and the calculated Cohen's κ coefficients are generally low and sometimes very negative (see Fig. 2 and Table S1). Due to the stronger signals of simulated winter MLD anomaly in the models, the RMSE reaches up to 668.71 m for GISS-E2-R and 639.89 m for KCM1.2.2, while the PMIP3 ensemble mean field has an RMSE of 112.51 m, comparable to that of the PMIP4 ensemble. Some models contributed to both the PMIP3 and PMIP4 MH experience with different model versions. Substantially improved agreement with proxy data from PMIP3 to PMIP4 is found for the GISS and the MRI models.

3.3 Transient simulations of the last deglaciation

The comparison of Holocene (the last 12 ka) time series shows little agreement among the models and between the models and the reconstruction (Fig. 6). The distribution of meltwater flux, either global and uniform as in HadCM3_uniform or routed to individual outlets as in HadCM3_routed, appears to have little influence on North Atlantic deep convection. TraCE-21ka remains nearly constant throughout the Holocene in deep-water formation regions and only shows a long-term trend in an area that is not a major deep-water formation site (the central subpolar North Atlantic). The MPI modelling group tested two ice-sheet reconstructions, Glac1D and Ice6G, with three model versions (P1, P2, P3). The results suggest that different ice-sheet evolutions generally lead to changes in timing and magnitude of MLD variations at the sub-millennial scale, without affecting the longer-term trends. In contrast, different tuning of the same model strongly affects the MLD response, as P1 differs substantially from the otherwise similar P2 and P3. This is expected given that P1 is an older version of MPI-ESM-CR and that P2 and P3 are equivalent to each other with a slight change in the threshold between cloud water and ice (Text S1 by Kapsch et al., 2022).

Figure 6 Centennial means of annual MLD at proxy sites from simulations of the last deglaciation, compared with proxy reconstructions of annual MLD averaged in 500-year bins. Proxy reconstructions are averaged in 500-year bins because many records are of low temporal resolution.

Cohen's κ coefficients calculated from the MH anomalies suggest that none of the models has a significant agreement with proxy reconstructions of annual MLD anomaly (Fig. 7 and Table S1). The RMSE relative to the reconstructions is very low for HadCM3_uniform, HadCM3_routed and TraCE-21ka, primarily because they simulate little variation over the last several millennia and hence very small anomalies.

Figure 7 Annual MLD anomalies (MH-PI) in the North Atlantic from the PMIP4 last deglaciation simulations (background) versus proxy data (dots).

4 Discussion

4.1 Inter-model differences

The inter-model spread of the North Atlantic MLD is significantly reduced across the PMIP4 ensemble compared to the PMIP3 ensemble (Fig. 5b, d). The PMIP3 ensemble mean field shows large anomalies and uncertainty in the central and eastern subpolar North Atlantic, which are not key areas of present-day deep-water formation. This suggests that the difference between the PMIP4 and PMIP3 ensembles is likely due to improved model physics, such as vertical mixing schemes and vertical grid resolution near the surface, which lead to more accurately simulated locations of deep convection in the new generation of Earth system models. In fact, improved mixing parameterisation is among the most important updates to the CMIP6-PMIP4 models (Xu et al., 2025; see also their Table 1). In this section, we will focus on differences among the PMIP4 models, which represent the current state of the art.

An important part of the changes in winter MLD simulated by the PMIP4 models seems to be tightly related to changes in winter sea ice concentration in MH relative to PI. Increased winter MLD is usually found at locations of reduced winter sea ice concentration (see Figs. 1 and S2), probably as a result of the two-way feedback between sea ice and open-ocean convection. A consistent reduction in Arctic sea ice minimum extent is found in the MH experiment across the PMIP4 models, primarily due to enhanced summer melt in response to higher insolation in Northern Hemisphere summer (Brierley et al., 2020). This leads to reduced export of Arctic sea ice in winter and reduced sea ice concentration in the subpolar North Atlantic, thus enhancing the ocean buoyancy flux to the atmosphere and promoting vertical mixing, which deepens the MLD. A more vigorous vertical mixing also allows warm and salty subsurface water of subtropical Atlantic origin to be mixed upward, which counteracts sea surface cooling and limits local sea ice formation, hence further strengthening the buoyancy flux and convection (Gelderloos et al., 2012).

Most models simulate a decrease in winter sea ice in different parts of the subpolar North Atlantic. In the Nordic Seas, this generally drives convection toward higher latitudes and the new sea ice edge (when sea ice extent exhibits a clear change). In the Labrador and Irminger Seas, the inter-model spread is large as the area ranges from entirely ice-covered to ice-free in the PI control experiment (Fig. S3). Consequently, when forced by the MH boundary conditions, the models produce remarkably different responses in the sea ice cover of this area. The initial state of sea ice concentration in PI can affect how its changes relate to the MLD. Unlike the previously mentioned negative correlation between MLD and sea ice concentration, when an increase in sea ice takes place in already densely sea ice-covered areas, there is no evident change in MLD since deep convection was already inhibited in PI and MLD remains very shallow (e.g., EC-Earth3, EC-Earth3-LR, INM-CM4.8, MIROC-ES2L). In contrast, in areas that are ice-free in PI, the MLD generally increases in response to enhanced ocean heat loss under reduced winter insolation, with no accompanying change in sea ice (e.g., CESM2, GISS-E2.1G, MPI-ESM1.2-LR, MRI-ESM2.0, NorESM2-LM). The intensified deep convection reduces stratification, so conditions there remain unfavourable for local sea ice formation.

Inter-model differences in sea ice and its transport seem to account for part of the large spread in MLD anomalies across the PMIP4 models. The models not only simulate considerably different pre-industrial sea ice extent as a baseline, but they also exhibit divergent responses to the MH forcings. The analysis of North Atlantic deep-water formation in the CMIP5 historical simulations (same generation as the PMIP3) suggests that deep convection is most realistic in models with realistic sea ice extent (Heuzé, 2017). Uncertainties in the PI sea ice extent are likely a more important source of the inter-model spread in MLD than uncertainties in the response of sea ice concentration. Our statistical analyses on winter sea ice concentration (Table S1) show no universal correlation between model performance in simulating MLD anomaly and sea ice anomaly: significant agreement on sea ice concentration anomaly does not guarantee agreement on MLD anomaly (e.g., FGOALS-g3, MPI-ESM1.2-LR), although the model with the best performance for MLD anomaly (MRI-ESM2.0) does have significant agreement for sea ice anomaly. This highlights the importance of other mechanisms affecting simulated MLD.

The sea ice model itself is affected by other model components. In some cases, differences in simulated sea ice and MLD may result from other model components. For example, by comparing the three versions of the EC-Earth3 model (EC-Earth3, EC-Earth3-LR, EC-Earth3-Veg), we find that using a lower-resolution atmospheric grid can cause a change in the amplitude of MLD and sea ice anomalies but preserves the spatial pattern, whereas incorporating a dynamic vegetation module instead of using prescribed modern vegetation leads to a marked change in the Labrador Sea and the Southern Ocean. Dynamic vegetation feedbacks enhance warming in the Sahara under MH forcings and increase rainfall in the Sahel region (Berntell and Zhang, 2024). This increased precipitation over northern Africa would be compensated by decreased precipitation over the western subtropical North Atlantic, creating a high salinity anomaly that is advected toward the subpolar North Atlantic (e.g., Zhang et al., 2021). In the EC-Earth simulations, this high salinity anomaly appears to affect the Labrador Sea the most, reducing stratification and sea ice formation, promoting convection and deep-water formation, and possibly enhancing the AMOC, which further strengthens the feedback loop. Vegetation can also increase warming in the Arctic through albedo feedback, leading to amplified Arctic sea ice retreat and reduced sea ice export in winter (e.g., O'ishi and Abe-Ouchi, 2011), which reinforce the stronger wintertime convection in the Labrador Sea. As another example, IPSL-CM6A-LR shares the same ocean and sea ice models as the EC-Earth3 suite, with a similar grid resolution but coupled to different atmosphere and land models. The simulated changes in MLD and sea ice are in opposite directions across the Nordic Seas between EC-Earth3 and IPSL-CM6A-LR. Similarly, FGOALS-f3-L and FGOALS-g3 produce considerably different MLD change and MH sea ice extent despite having the same ocean and sea ice models, implying the governing role of the atmosphere component in a coupled model as the primary driver of open-ocean convection and sea ice formation.

These findings highlight the importance of air-sea-ice interactions in modulating convection in the subpolar North Atlantic and the need to resolve both oceanic and atmospheric processes, as well as sea ice dynamics, to correctly represent the ocean mixed layer. In addition, the difference between HadGEM2-CC and HadGEM2-ES suggests that the implementation of tropospheric chemistry might mitigate the model's response to MH orbital forcing (see Fig. 2). However, with only 35 years available for HadGEM2-CC, we cannot establish whether this difference exceeds internal variability, and we note this as a potential direction for future work. Models differ substantially in the North Atlantic deep-water formation areas (Fig. 5b), indicating that accurately simulating the strength of deep convection and the MLD remains challenging. The situation is similar for simulated AMOC strength, with a large spread across the PMIP4 models (Brierley et al., 2020; Jiang et al., 2023). Although many models simulate little change in AMOC strength, the direction of change is generally consistent with the relative proportion of increase and decrease in winter MLD. For example, EC-Earth3-LR simulates a stronger AMOC in the MH relative to PI while IPSL-CM6A-LR simulates a slightly reduced one (Jiang et al., 2023).

4.2 Proxy-model discrepancies and influence of meltwater forcing

While most PMIP3 and PMIP4 models agree that winter convection in the MH is stronger than PI in the Labrador Sea because of the deeper MLDs, proxy data suggest the opposite (Fig. 5a, c). Differences between simulated and reconstructed sea ice change might explain the discrepancy at the proxy record site in the west and the other two sites in the south, but not in the northern Labrador Sea (Fig. 5e). In the northern Labrador Sea, both proxy records and numerical simulations indicate reduced winter sea ice concentration relative to the PI; the reconstructed and simulated sea ice concentrations are both below 50 % in the MH (Fig. S4). However, while models simulate increased MLD straightforwardly in this case, reconstructed MLD decreases.

This leads to the hypothesis that during the Holocene, meltwater influx played a key role in modulating deep convection and deep-water formation in the Labrador Sea. The MH experiment is an equilibrium simulation designed to examine the impact of changes in orbital forcing at times when greenhouse gas levels were similar to those of the preindustrial period, and the continental configurations were close to modern ones (Otto-Bliesner et al., 2017). Hence, the modern configuration of ice sheets is prescribed in the MH simulation, and the ice sheets are in equilibrium with other components of the Earth system in the final outputs, which is a known source of uncertainty (Otto-Bliesner et al., 2017). In the real world, the MH Labrador Sea likely remained under a superimposed effect of rapid ice-sheet melting before the Middle Holocene and ongoing ice-sheet melting, thus far from an equilibrium state. Geomorphological evidence indicates a final disappearance of residual ice caps and glaciers in eastern North America as late as 6–5.5 ka BP (e.g., Clark et al., 2000; Dyke, 2004; Dalton et al., 2023), and glaciological, ice cores, as well as paleoecological records from coastal areas agree that the Greenland Ice Sheet reached a minimum extent around 4000 years ago (Briner et al., 2016 and references therein, including Vinther et al., 2009). Marine proxy records indicate highly variable sea surface salinity in the Labrador Sea during the Early to Middle Holocene (e.g., de Vernal and Hillaire-Marcel, 2006; Gibb et al., 2014), with lower sea-surface salinity until the Late Holocene at some locations (e.g., Allan et al., 2021; Solignac et al., 2004), which could be attributed to meltwater discharge. Proxy data of surface and subsurface waters in the Labrador Sea suggest the presence of a buoyant surface layer unfavourable for deep-water formation until the Middle Holocene (Hillaire-Marcel et al., 2001; de Vernal and Hillaire-Marcel, 2006), suggesting possible changes in the routing and rate of meltwater input. Together, the proxy evidence indicates that the Labrador Sea was under prolonged influence of meltwater influx throughout the Middle Holocene, which likely inhibited deep convection. Such transient effects and the continued meltwater flux into the North Atlantic are missing in the MH equilibrium simulation.

In contrast to the Labrador Sea, there is relatively good proxy-model agreement in the Nordic Seas despite the missing meltwater forcing, suggesting that the northwestern North Atlantic may be the most affected deep-water formation region under a future scenario of global warming and ice-sheet melting. This is consistent with modelling studies of meltwater pathways in the present-day North Atlantic that suggest generally limited exchange from the East Greenland Current into the Nordic or Irminger Seas, with most of the meltwater ending up in the Labrador Sea or Baffin Bay (e.g., Dukhovskoy et al., 2019; Pennelly et al., 2019).

In the last deglaciation simulations, models generally do not overestimate MLD in the Labrador Sea as much as in the MH equilibrium simulations (Fig. 6), which could provide evidence for the influence of meltwater fluxes. However, the models do not show better overall agreement with the proxy records or with each other, raising questions about the robustness of this evidence. The models' sensitivity to meltwater forcing appears highly variable, as simulations run with the same ice sheet reconstruction (e.g., HadCM3_routed versus the MPI_Ice6G ensemble) yield very different evolutions of MLD, with some remaining nearly constant and others showing large fluctuations. While reconstructions in the Labrador Sea reveal an increase of MLD toward modern-like conditions over the Late Holocene, none of the models reproduces this trend. This implies that the memory of the ocean in long transient simulations might be too long compared to what proxy records suggest, especially in the Labrador Sea. In the Nordic Seas, the MPI ensemble largely overestimate the MLD, which is consistent with the previously mentioned fact that meltwater affects the Labrador Sea more than the Nordic Seas; the HadCM3 simulations are the closest to the proxy records in the Nordic Seas; TraCE-21ka shows lower annual MLD than proxy records, and the results further indicate that this model is probably not resolving the correct locations of deep convection. The divergence between transient simulations and proxy data in North Atlantic deep convection could be one of the reasons behind the `Holocene temperature conundrum' discussed by Liu et al. (2014).

Unfortunately, the robustness of the statistical analyses in this study is limited by the small sample size. The estimated statistical metrics, notably the Cohen's κ coefficient, may not accurately reflect the true performance of the models and should therefore be interpreted with caution. Nevertheless, it is interesting to note that the ensemble average does not always align better with proxy data than individual models, and that models capable of, or nearly capable of, reproducing the direction of MLD change do not necessarily have a smaller RMSE. The response of MLD to MH forcing could be dampened in some models due to the apparent deep bias in some deep-water formation regions (Figs. 3 and 4). The model biases could be due to the various parameterisations, such as those for mesoscale eddies or vertical mixing scheme. The coarse resolution of the models could also introduce bias, as low-resolution models tend to exhibit larger deep biases in winter MLD in some eddy-rich regions than their high-resolution counterparts (Treguier et al., 2023). For transient simulations of the last deglaciation where a meltwater flux is applied, the coarse resolution models struggle to simulate the narrow, fast boundary currents, in addition to the model biases discussed above. The slower boundary advection would cause freshwater to be trapped within the ocean basins for much longer time. The biases could also be related to uncertainties in boundary conditions, such as the ice-sheet reconstruction and the associated timing, volume, and routing of meltwater, changes in vegetation cover, etc. For the PMIP MH simulations, using monthly rather than higher-frequency model output to recompute MLD would likely yield smoothed signals of density profiles and extreme events such as storms (Treguier et al., 2023), further adding to the uncertainty.

The model biases are challenging to analyse here as there is uncertainty in the reconstructed MLD data as well (see Wu et al., 2025). The main sources of proxy-related uncertainties are errors in sediment sample age estimation and errors in the reconstruction method. The errors in sediment core age modelling vary from core to core and depend on the number of ¹⁴C dating points in each core and the accuracy of the estimated ocean reservoir-age correction. This error can be as high as a few hundred years, affecting the actual 6.5–5.5 ka BP time window used to represent the MH in sediment cores. Its impact on the results, however, is difficult to analyse and core-specific. The reconstruction error (RMSEP) is estimated from leave-one-out cross-validation to be 40.9 m for winter MLD, 16.0 m for annual MLD, and 17.6 % for winter sea ice concentration. Finally, we found that the reconstruction method can exhibit regional biases, as discussed in Wu et al. (2025), which we attempt to mitigate using the MH minus PI anomaly approach. For example, the reconstruction method tends to underestimate the MLD in the southern Labrador Sea and the Iceland Basin (see Fig. S12 in Wu et al., 2025). Overall, the proxy-related uncertainties are probably comparable to, or slightly less than, the uncertainty of the multi-model mean. Systematic benchmarking of the previous generation of PMIP MH experiment shows that differences in model performance are only weakly related to modern-day biases (Harrison et al., 2015). Greater availability of quantitative reconstruction data based on different proxies is therefore necessary to improve climate model evaluation.

5 Conclusions

Our results show that the models generally simulate a realistic seasonal cycle of the MLD in key deep-water formation regions, in the sense that the MLD slowly deepens in late autumn, reaches a maximum in late winter, and then rapidly shoals in spring. Compared to the PMIP3 ensemble, the PMIP4 ensemble shows substantial improvement with reduced inter-model spread and more realistic locations of deep convection in the North Atlantic (Fig. 5). However, the simulated response of winter MLD to MH forcing still varies significantly across the PMIP4 models, and the proxy-model agreement appears to be region-dependent. While there is relatively good agreement in the Nordic Seas, significant discrepancies in the Labrador Sea persist from PMIP3 to PMIP4. This confirms that the discrepancies in the Labrador Sea are not artefacts of a single model generation, but rather a systematic bias in the MH simulations. We propose that the missing meltwater forcing in the PMIP midHolocene simulation accounts for this marked discrepancy between proxy data and model simulations, implying a higher sensitivity of the Labrador Sea to changes in meltwater fluxes in the subpolar North Atlantic. Comparison against simulations of the last deglaciation lends some support to this hypothesis, in that the overestimation of Labrador Sea MLD is weaker in these experiments (Figs. 3, 4, 6). However, confidence in this interpretation remains limited by uncertainties in the spatial and temporal distribution of the applied meltwater flux, and by biases in the deglaciation simulations themselves. The difficulty of discerning the impact of meltwater during the MH therefore reflects the uncertainties of current climate models and experimental designs (Sect. 4.2), rather than constituting evidence against a role for meltwater.

Despite inherent uncertainties in the proxy data and potential model biases, the present study on the MLD demonstrates the need to improve the representation of key processes and feedback mechanisms at the atmosphere-ocean-sea ice interface, which play a determinant role in driving open-ocean convection. This is particularly critical for the subpolar North Atlantic, a potential tipping point area (Brovkin et al., 2021). Further investigation will require Holocene experiments with improved meltwater forcing in the latest generation of PMIP-CMIP models.