This is a good paper that presents some interesting new simulations. I appreciate the work that's gone into these runs and their analysis and can readily see this manuscript being published in Climate of the Past. There are some aspects of it that need clarification before publication, and I think a little bit of further analysis would greatly enhance the reach of this manuscript. I especially appreciate the data and code placed in the online repository.

1. The model description (Sect 2.1) mentions nothing about the land surface model. Given the importance of the vegetation fraction in this manuscript, you need to provide some information about how vegetation is simulated by the model (tree, grass etc) – and what, if any, feedbacks it has on the atmosphere.
2. I feel the analysis about the rates of change (Sect 3.3) is out of place in this manuscript. It seems to invoke a fundamentally different conception of an AHP to the other work. The rest of the work talks about thresholds (implying transitions between bistable states). Yet this section discusses the speed of the changes as being related to the speed of forcing changes irrespective of their location w.r.t. the thresholds. Personally, I feel this aspect of the research should be removed to focus more on the subject in the title.
3. You discuss the threshold as a function of the maximum orbital forcing. This may be appropriate for precipitation, but is this really the best way to think of vegetation threshold? Intuitively, I see a threshold as being lower than the maximum value with the intensity of the vegetation response driven by the time spent over that threshold.
4. It is not clear precisely what is plotted in the trajectories of simulated data. Are these the data for a single grid box? If so, which one? Is the vegetation fraction presented a proportion of this grid box, with the rest of it being bare soil?
5. Why have you selected only the past 190 kyr (Sect 2.2)? I presumed this was motivated by the 2 references cited on L36 – although you should make this explicit. It seems though that Ehrmann & Schmiedl review back to 200ka and Blanchet et al seems to go back to 160ka from their Fig 3. I don’t expect you to redo any simulations – your start date is fine for the science. But it needs a solid motivation written in the paper.
6. There is no discussion in the paper of internal variability in the simulations. My own work (Brierley et al, 2018, https://www.nature.com/articles/s41467-018-06321-y) building of Zhengyu Liu’s model relies quite heavily on the fact that the AHP transitions involved some stochasticity. I suspect this will be case for CLIMBER-2 as well, and that this would explain the difference in precipitation at MIS5e between EI7 and E0 in Fig 6b. Again, I don’t think any additional analysis is needed – just some discussion of its implication for your analyses.
7. You could go further with your simulations and combine the results from the future simulations with that of EI2, EI4 and EI6 to perform an analysis similar to that in Fig. 3 to quantify the impact of GHG forcing on the orbital threshold. As currently written this feels like a missed opportunity to really demonstrate the statement in the title.

Other comments:

• ‘Synergical’ feels very awkward – try ‘synergistic’
• I agree that with Dr Liu that a slight rebranding of the Monsoon Index would be helpful
• You should explain how the lagged peaks in Fig 2a reflect the intensity during the sapropel. You make no comment about the split event at 5c in SL77. Why are these better measures of intensity than something like the co-eval Ba/Al ratios measured by Zeigler et al (2010)?
• “reckon” on L169 sounds informal. Please replace.
• You are too precise stating that the change point at 20Wm-2. Surely all you can tell is that its between 15-20 Wm-2.
• How do you justify LOWESS smoothing all the forcing in Fig. 4, but not the simulated vegetation fraction? [I recommended cutting this section above]
• I strongly suspect that the analysis in Fig 5 would have also show the rates of initiation and termination of the AHP events is strongly correlated to the peak monsoon index. How can be sure that your style of analysis is more appropriate. [I recommended cutting this section above]
• 6. I like this figure, but can you please check that it works for color-blind individuals.
• This sentence seems odd. If you really feel that it is only the weak orbit that matters, then please rephrase to avoid the conflation with ‘glacial times’ – as that phrasing intuitively suggest that GHG and ice-sheets play a role. You might want to try: “This analysis demonstrates that it is the relatively low maximums in orbital forcing that result in the absence of AHP conditions at 6b, 4 and 3a – rather than the low GHG forcing or large ice sheets.”
• It would be instructive to take the work about future AHP conditions a little further. Can you find a way to quantify the impact of GHG forcing on the orbital threshold. I feel that there should be enough data here.
• I also wonder if you could provide some additional context for the future simulations for those of us not fully versed with the future carbon cycle pulses. As well as the GHG forcing, it might be helpful to plot global mean temperatures and atmospheric CO2 levels. In effect, I am wondering how the future AHP at M1 relates to proposed warming levels and safe operating spaces.

This study presents modeling results of the last 190,000 years of African rainfall and vegetation history, classifying certain thresholds as “African humid periods” and commenting on the strength, duration, and rates of change of these differing AHPs. The authors find that orbital forcing is the primary driver for changes in rainfall and vegetation extent during past AHPs, but that the sensitivity threshold of AHPs to orbital forcing is modulated by GHG concentrations. Future modeling experiments are also conducted that show future AHPs are more likely to occur with higher concentrations of GHGs, as future orbital insolation thresholds are too low to induce AHPs without GHG increases.

This study is well motivated and provides novel findings with respect to previously unknown factors contributing to the strength, duration, and rates of change of past AHP. This paper is exceptionally well-written and clearly presents its results and conclusions. In addition to a few minor comments, I believe one area for improvement can come from some added discussion on the uncertainties present within the very coarse model resolution of CLIMBER-2. It is important to show that the authors have considered all of the uncertainties involved in using this specific model and conclude that these uncertainties do not impact the conclusions of this paper – i.e., this model is the perfect fit for use with this specific research question.

I recommend this paper be accepted with some very minor revisions. I list each comment for the revised manuscript below.

Comments:

-It will strengthen the manuscript to elaborate upon the scale of the research question with regard to these simulations (for example: this study examines shifts in state of climate, such as desert vs. >50% vegetation cover, present within the single North Africa grid cell and does not require finer details with regard to the simulated climate) and how examination at this scale minimizes the large uncertainties present with using CLIMBER-2 to simulate paleoclimate. Bringing in discussion of multiple climate equilibria (green vs. desert) in northern Africa may help to strengthen this argument.

- In Table 1, it would be more clear to list “Monsoon index via orbital parameters” (or something like this) so to not confuse readers over what is being prescribed in the model. The authors prescribe orbital parameters, which in turn dictate the monsoon index, rather than directly prescribing “monsoon index” as a specific boundary condition. Slight added nuance to reflect this would preclude confusion for future readers.

-In Table 1, what does GHG radiative forcing = 0.0 W/m2 correspond to? The base value is listed for monsoon index (line 428), so it would be helpful to include the same for GHG radiative forcing. Or if this value is more difficult to assess, at least define more clearly that deltaRF is a change from the modern day… which is what time period? 1950 CE?

-On line 158, there are several studies with updated simulations using sophisticated models that could be cited here, in addition to Harrison et al. (2015). I would suggest adding at least a few of the following citations:

Pausata et al. (2016), 10.1016/j.epsl.2015.11.049

Thompson et al. (2019), 10.1029/2018GL081225

Hopcroft et al. (2021), 10.1073/pnas.2108783118

Chandan & Peltier (2020), 10.1029/2020GL088728

Dallmeyer et al. (2020), 10.5194/cp-16-117-2020

-Both the interglacial and glacial factor separation analyses are important results of this paper, yet only one is presented in the main text. I would suggest the authors bring the glacial factor separation analysis into the main text as an additional figure. Or the authors could at least describe why they believe the interglacial case is more important than the glacial case and use this explanation to justify why the interglacial case is included in the main text while the glacial case is not.

This paper is a valuable contribution to the theory underlying African Humid Periods and their variable forcings. The authors present a carefully considered set of intermediate complexity model simulations that allow for factor separation analysis. They have clearly produced a lot of data/results, and I appreciate the efforts they have made to condense the work to the most important points. I think the paper is close to being ready for publication. Here, I touch on some previous points by other reviewers that I agree with, and I add a couple of additional, minor points.

The most substantial point that I wish to emphasize comes from Dr. Liu about how the main text is somewhat disconnected from the title. This is a substantial point only in that I think the paper could benefit from re-structuring the arguments, but I don’t see this as necessary for publication. Specifically, I suggest re-framing the paper more explicitly as a comparative analysis of past and future AHPs. This would involve discussing the future simulations more prominently and, as other reviewers mentioned, diving into more detailed hypotheses as to how/why GHGs lower the orbital threshold. The question of whether emissions can compensate for low future eccentricity to cause future AHPs is thought-provoking, and the results—casting emissions scenarios as the primary determinant of the frequency and amplitude of future AHPs—could motivate much further research into GHG and orbital “synergies”.

A couple of smaller points that I agree with from other reviewers. I like Dr. Liu’s suggestion to call “Monsoon Index” the “Monsoon Forcing Index”. I also agree with Dr. Brierly that more background on the land surface model is needed, especially with respect to the threshold behavior, relevant feedbacks (including fire), and whether there is hysteresis. I also agree that the rate of change analysis could be removed. It’s not currently grounded by anything in the discussion, and I agree that it is difficult to square with the threshold behavior.

My two suggestions are (1) to cite/discuss some more proxy work; and (2) be more explicit about any assumptions associated with factor separation analysis:

1. The paper focuses on two records from the Mediterranean for comparison. However, there are other records that span the same time interval, and it is worth mentioning how they compare (amplitude, duration, etc) to the new model results. Two datasets that are particularly relevant are Miller et al. 2016 (JQS) and Skonieczny et al. 2019 (Sci. Adv). The Miller paper is useful for more directly comparing the vegetation results to a vegetation reconstruction; and the Skonieczny paper presents another dust flux record off West Africa.
2. One concern I have has to do with any assumptions inherent to the FSA (I don’t have expertise in FSA, so please bear with me). It seems like one implicit assumption is that any non-linearities (when multiple forcings yield a different result than the sum of individual forcings) can only arise due to “synergies” or interactions between the forcings. That is, the response to any forcing is assumed linear so the responses can be summed together (and deviations from the sum are synergies). However, a non-linear response to a forcing (such as threshold behavior in vegetation %) could lead to “apparent synergies” between forcings that are actually projections of a non-linear response (rather than a non-linear interaction between forcings). For example, if GHGs and orbital forcing alone both cause a stepwise increase from desert-to-grassland, then FSA would expect the GHG + orbital simulation to produce this stepwise transition twice (without synergies) and any deviation from this would be attributable to synergies. The authors briefly touch on this specific case in line 265, but they relate the issue to the fact that vegetation % is bounded between zero and 100. However, maybe it’s the case that the bounding only makes the issue clearly diagnosable. I’m curious if the issue is broader, applying to any non-linear response where the response is not necessarily the sum of its parts (even without “synergies”). I expect that a basic discussion of the assumptions of FSA would suffice here. If my comment about non-linear responses is entirely off-base, then maybe adding a sentence about why this intuition is wrong could be helpful for readers like me.

References:

Miller, C. S., Gosling, W. D., Kemp, D. B., Coe, A. L., & Gilmour, I. (2016). Drivers of ecosystem and climate change in tropical West Africa over the past∼ 540 000 years. Journal of Quaternary Science, 31(7), 671-677.

Skonieczny, C., McGee, D., Winckler, G., Bory, A., Bradtmiller, L. I., Kinsley, C. W., ... & Malaizé, B. (2019). Monsoon-driven Saharan dust variability over the past 240,000 years. Science advances, 5(1), eaav1887.