Vaittinada Ayar et al. examined a divergent projection of the ENSO-CO2 flux relationship between CMIP6 ESMs in the SSP5-8.5 scenario. Half of the ESMs simulate a reversal ENSO-CO2 flux relationship while the other half EMSs have a consistent ENSO-CO2 flux relationship from historical to future period. They found the reversal relationship in the first half ESMs is due to a faster-increasing surface DIC, enhanced primary production variability, and pCO2 increase. The manuscript is well written and well structured. These findings are very interesting and have significant implications for understanding the different CMIP6 model projections and model bias. However, I have some concerns about the missing terms in the analysis of the contribution to the CO2 flux and ocean pCO2 variability. I would recommend it for publication in ESD after the concerns and comments are considered by the authors.

General comments

1. Since the air-sea CO2 flux is related to three terms: ocean pCO2, air pCO2, and wind-solubility coefficient, the authors only analyze the ocean pCO2 in the manuscript. I am very interested in the role of wind-solubility coefficient and air pCO2 in explaining the divergence in two groups of ESMs. Different models might have different wind and temperature variability which might contribute to the CO2 flux variability. It is worth quantifying and discussing the wind and solubility terms. All the ESMs might use the same air pCO2, so the air pCO2 might have a very little contribution. However, it is needed to be at least discussed in the manuscript.

2. Ocean pCO2 is sensitive to four terms: temperature, DIC, alkalinity, and salinity (Takahashi et al., 1993). The authors only discuss the temperature and DIC. Although the DIC is dominant in the ocean pCO2 variability, the alkalinity has a very large compensation. The alkalinity might partly contribute to the model divergence. In addition, the precipitation probably changes a lot under global warming (Cai et al., 2015), this might drive a relatively large variability of alkalinity and salinity in the future. It would be convincing if the author could discuss/quantify the contribution of alkalinity and salinity to the model divergence?

3. Line 245-247. The authors found two differences between two group of ESMs (Large increase of surface DIC and lower range of DIC changes). Fig.9 could show the large increase in surface DIC. However, I could not see a figure showing the lower range of DIC changes during ENSO phases. I would suggest one such figure in the main text or supporting information.

4. Fig. 9. Except for surface DIC difference between preserved and reversed models, I also see the difference in subsurface DIC (e.g., 200-300 m) between two groups of ESMs. What is the role of subsurface DIC difference in the model CO2 flux-ENSO relationship divergence? Why is the subsurface DIC also different in the two groups of models?

Minor comments:

1. Line 41-43. The tropical Pacific ocean CO2 flux anomaly is not only related to the upwelling strength but also related to the poleward Ekman transport driven by easterly trade wind. One more reference (Liao et al., 2020 GBC) is suggested.
2. Line 80. What is the re-grid method?
3. Line 91. The ENSO variability is usually an interannual variability ranging from 3 to 7 years. Could the author plot the total CO2 flux and CO2 flux anomaly at a sample point to show how well is the detrend method? Could the detrend method remove the decadal variability?
4. Line 100. Why do you use 1981-2010 as the climatological period instead of 1985-2014 which is the contemporary period defined in the manuscript.
5. 2 Caption. What is the observed data of SST and CO2 flux? I know the authors state them in the method section. However, it would be clearer for the readers if the authors could detail them in the caption. For example, SST is JRA.
6. 3. Why do the authors use 5N-5S instead of 2N-2S? This is not consistent with the method section.
7. What is the CMIP6 ensemble anomalies one standard deviation?
8. Line 379. The text reads like Liao et al. (2021) selected the model subjectively and got a partial and biased conclusion. Actually, Liao et al., (2021) use a strict and reasonable constrain method to select the model. The results are physically rational and convincing. I would suggest rephrasing the words. A suggested way would be: “With a strict constrain method based on contemporary observations, the model tends to show a weaker future CO2 flux anomalies during ENSO phases (Liao et al., 2021).”
9. Lines 381-382. The increasing Revelle factor and ocean pCO2 sensitivity to temperature would be a general result in my opinion. This point is discussed by many studies. I would rephrase or delete this point.

References:

Cai, W., Borlace, S., Lengaigne, M., van Rensch, P., Collins, M., Vecchi, G., et al. (2014), Increasing frequency of extreme El Niño events due to greenhouse warming, Nature Climate Change, 4, 111.

Liao, E., Resplandy, L., Liu, J., and Bowman, K. W. (2020), Amplification of the Ocean Carbon Sink During El Niños: Role of Poleward Ekman Transport and Influence on Atmospheric CO2, Global Biogeochemical Cycles, 34(9), e2020GB006574.

Takahashi, T., Olafsson, J., Goddard, J. G., Chipman, D. W., and Sutherland, S. C. (1993), Seasonal variation of CO2 and nutrients in the high-latitude surface oceans: A comparative study, Global Biogeochemical Cycles, 7(4), 843-878.

Firstly, Vaittinada Ayar et al. evaluated the skill of 16 ESMs from CMIP6 to reproduce the ENSO-CO2 flux relationship in the Pacific Equatorial. For this, model outputs of chemical, physical and biological variables were compared to observational datasets. Secondly, they analysed how the simulated ENSO-CO2 flux relationship evolved in the future (i.e., 2071-2100) using model projections under the SSP5-8.5 scenario. They found that half of the ESMs projected a positive correlation between ENSO-associated warming and sea-air CO2 flux anomalies, instead of the current negative correlation. According to their findings, the future reversal of the ENSO-CO2 flux relationship is induced by the thermal component of pCO2 becoming more important than the non-thermal components. However, they concluded that this reversal is unlikely because it could be related to model biases in the historical period.

Such comparative study of CMIP6 models improve our understanding of the influence of climate modes on the carbon cycle and could provide useful metrics to evaluate ESMs. The manuscript is clear and well written. However, I have some questions about the evolution of the simulated CO2 flux anomaly variability which is different between ESMs identified as “reversed” or “preserved”. Therefore, the paper will likely be a significant scientific contribution with minor revisions.

General comments:

1) Figure 5 and Line 206: “This reversal is thus independent on the performance of the model over the contemporary period, though the models in the first row tend to simulate lower than observed CO2 flux anomaly variability.”

Firstly, after reading the manuscript, the results suggested that the reversal behavior was indeed induced by the model performance in the contemporary period. Authors should modify or clarify this sentence.

Secondly, when looking at figure 5, the lower CO2 flux variability in the “reversed” ESMs than in the “preserved” ESMs is a striking feature. I would like to see some discussion about the influence (or the relationship) of this feature with the conclusions. For example, could the historical low CO2 flux variability in the “reversed” ESMs be related to their higher carbon uptake than in the “preserved” ESMs? Authors focused on the understanding of the correlation between the annual CO2 flux and the ENSO index, but could some of their findings explain the variability in the amplitude of the simulated CO2 flux anomalies? As a reminder, most models underestimated the CO2 flux variability (line 197 and Table 3) and according to the figure 5 this is more visible in the “reversed” ESMs.

2) Authors estimated the depth of the thermocline (line 105) but their discussion and conclusions focused on the stratification (or the vertical gradient), which are two different concepts. Although there is no difference between the two ESM groups in term of “thermocline depth” (line 306) the vertical stratification might be different. Therefore, could authors replace their “thermocline depth” estimate with a stratification estimate.

Minor comments:

3) Line 28: “…the Equatorial Pacific CO2 flux represents the dominant mode of variability of the global oceanic CO2 flux variations (Wetzel et al., 2005; Resplandy et al., 2015…”. According to Resplandy et al. (2015), for some ESMs, the Southern Ocean can also be the dominant mode of variability of the global oceanic CO2 flux variations.

4) Line 110 – At which temporal resolution is the thermocline depth estimated? Monthly?

5) Line 175: “Note that the observed average is the result of the climatology over the 2004-2017 period while the average for CMIP6 is computed over 30 years (1985-2014).” Could authors calculate the CMIP6 climatology using the same period (i.e., 2004-2017)? If not, this information should be included in Figure 3.

6) Line 188: “The correlation between annual CO2 flux anomaly and annual ENSO index is given for the models for each 30-year sliding window over the 1850-2100 period.” Why did author choose a 30-year sliding window? Is it the observational period? This information needs to be added.

The study by Ayar et al. aims to investigate the ENSO-induced response of air-sea CO₂ fluxes from CMIP6 ESMs. ENSO is a big driver of inter-annual variability in the carbon cycle, and studying the response of the carbon cycle to ENSO is critical for understanding how feedbacks between the physical climate system and global carbon cycle operate. The study explores the response of our current generation of ESMs under a high CO₂ future climate scenario (future period 2071-2100) in the context of their skills over a contemporary reference period (1985-2014). The authors use a wide variety of observational data and skill metrics to show distinct relationships emerging from the different ESMs. The study is timely and significant and should be considered for publication once the authors address and clarify a few issues for the benefit of the reader and the scientific community. Specifically, I would strongly advise the authors to consider and respond to Major Comment #1 given that the authors have submitted to ‘Earth System Dynamics’ - a journal that does emphasize on interdisciplinary Earth system science, beyond just ocean sciences.

Major Comments:

1. Even though the authors are looking at the response of coupled models, the authors have ignored any changes or even providing  statements about the atmospheric response (except Lines 340-344). There is no mention of changes in trade winds and/or changes in conditions of the air-sea interface due to the weakening of the easterly trade winds (during El Niño, for example). To me, this is a key ingredient that is missing from the study. This is a CMIP-based study and since the tropical ocean-atmosphere are strongly coupled with each other, the authors do need to provide qualitative statements about how atmospheric conditions across the ESMs (preserved vs. reserved) evolve that impact the oceanic ENSO response. Quantitative analyses regarding changes in atmospheric winds across the study time periods (or a figure or two) would be better, but I recognize that a quantitative evaluation of dynamical wind response is not a trivial task.
2. The authors should consider evaluation of the models for specific ENSO cases - strong El Niño or strong La Niña years. Figure 5 provides a first indication that the “preserved” ESMs agree better with the observations than the reversed ESMs. But the comparison is noisy, and it may be better to examine specific strong and very strong ENSO events between 1950 – 2014. Approximately 10 such events can be identified for both El Niño and La Niña conditions that should allow robust statements on which of the two groups of ESMs (preserved vs reserved) validate better against observations.

Minor Comments:

1. Line 2 – change to ‘over the tropical ocean less carbon is released during El Niño…’
2. Line 56 – it is hard to interpret what the authors mean by the phrase ‘an end member future projection’. While it becomes clear eventually that the authors are referring to the high-emission scenario, maybe add a sentence or two here to clarify this phrase for the benefit of the reader.
3. Line 97 – please check the grammar and punctuation
4. Lines 145-146 and Lines 337-339 – it is a bit strange that while the authors define a classical Niño 3.4 domain (Lines 99-100), the study area is subsequently shifted to different longitudes. This matters because not all El Niños are similar and whether we are looking at an EP or a CP El Niño should have implications for the findings of this study. Did the authors consider evaluating the model simulations based on different El Niño types?
5. I would strongly encourage a modified version of Figure 11 – again, instead of looking at 1850-2100, maybe pick a period or specific strong & very strong ENSO years, for which the authors can plot a ‘best estimate’ of air-sea CO₂ flux from observations and/or models (for example, see Ishii et al., 2014, Biogeosciences, https://doi.org/10.5194/bg-11-709-2014). It would be interesting to see which of the ESMs actually fall in the region where surface CO₃ ²⁻ concentration obs. and the most optimal estimate of air-sea CO2 fluxes overlap. Can we identify a subset within the 16 CMIP6 ESMs that validate better against the observations? This study has already laid the foundation for providing this key message, thereby really helping the improvement of future ESMs and CMIP simulations.