This manuscript explores methodological aspects of rare event simulation, which is a promising line of research to improve our understanding of extreme events in a broad range of complex systems.

Specifically, the authors study in a systematic manner some important design choices for resampling algorithms applied to deterministic systems such as the timing and structure of the imposed perturbation and their impact on the ability of the algorithm to sample more extreme events.

According to me, the major novelties in the manuscript are the following:

- introducing an estimator for the probability of resampled trajectories in the "ensemble boosting" method, which allows to compute the climatological probability of the rare events under study.

- introducing a different method to generate perturbations at resampling times compared to what was done before in the litterature, constructed by sampling a low-dimensional space. This allows to study in details the relation between the perturbation and the resampled amplitude of the event (called severity in the manuscript).

- introducing optimality criteria to select "advance split times".

To investigate these methodological aspects, the manuscript uses as an example the fluctuations of a tracer concentration in a baroclinically unstable two-layer QG model.

The work presented in this manuscript is very rigorous and it is presented in a precise manner. The problems under study are an important part of methodological developments which have a great potential for wide applications. For this reason, the manuscript should be very useful for researchers wishing to implement rare event algorithms in practice. In my opinion this is very high-quality work.

Perhaps the only general shortcoming I can see is that due to its technical nature, the manuscript is a bit difficult to read, and its conclusions are restricted to methodological aspects. Maybe one way to broaden its potential audience could be to reinforce the physical content, for instance by discussing the potential importance of extreme mixing events in baroclinic turbulence (typical theories diagnose effective diffusivity based only on typical fluctuations). Given the impressive amount of work that already went into the current manuscript I would not make this a condition for recommending acceptance of the manuscript, but only a suggestion that the authors might want to consider to extend the geophysical impact of their work.

In addition to this general comment I have a few questions on the technical content of the manuscript.

- The question of the effect of introducing a random perturbation on the statistics is very interesting. In the example studied in the paper, Fig. 13 shows that some AST selection procedures leads to apparent bias. Presumably this is just a sampling issue: because these AST selection procedures are less efficient at generating larger severities, the tail probability tends to be underestimated. In addition to the empirical evidence, can you obtain analytical insight on the bias of your estimator?

- Could you extend the method to estimate statistics of any observable conditioned on the fact that the severity is above a given threshold within ensemble boosting? In diffusion Monte-Carlo algorithms this is possible because it performs importance sampling in trajectory space. Is it the case here?

- another interesting aspect of the manuscript is the idea to use construct perturbations from a low-dimensional sample space. Intuitively it seems that it should lead to smaller variance of the estimators. However, in the case considered here I wonder if this would necessarily be the case. If the growth of the perturbation is indeed governed by the linear baroclinic instability mechanism, a random perturbation directly constructed in streamfunction or vorticity space should project onto the most unstable mode (the same one you are forcing by construction in the manuscript), and the evolution should quickly be dominated by this mode. Could you comment on this?

- while the optimality criteria introduced in the manuscript are interesting, for practical use one would have to estimate the AST online, without systematically searching for the optimum. Do you expect this to be a limitation for applications?

- I understand that it is not the direct goal of the paper to design an algorithm which is more computationally efficient than DNS. Nevertheless, it would be natural to expect that it should be the case even without any specific effort. In previous uses of rare event algorithm with climate models for instance, there was an immediate gain of orders of magnitude. Here, Fig. 13 suggests that ensemble boosting with proper AST selection indeed performs better than equal-cost DNS, but it is difficult to appreciate how much really. Would it be possible to do such a comparison, for instance in terms of the return time of the most extreme events simulated? Unlike algorithms of the interacting particle systems family, ensemble boosting is not iterative: the initial ancestors are used to try to simulate more extreme descendents, but these descendents themselves are never used as ancestors again. Do you think this plays a part for computational efficiency?

Specific comments:

- some of the questions formulated in the abstract are very general and the manuscript addresses only part of them. I would be in favor of a more focused abstract.

- Fig. 1: the caption is long and it is not immediately evident when looking at the figure what is plotted on panels b)c)d)iii; perhaps add a legend directly on the figure to make it more clear? Even with the caption it is not completely transparent why there are two solid red lines in those panels.

- p10: "can be estimated it by...": remove "it". 

- p18, paragraph 2: the relevant timescales, in particular the eddy turnover time, is estimated empirically here if I understand correctly. Is it compatible with the phenomenology of baroclinic turbulence, which should allow you to estimate it a priori from model parameters?

- p18, paragraph 3: the results suggest that mixing is more efficient in eastward jets, can you relate this to known results on the phenomenology of baroclinic transport?

- Fig. 4: shouldn't the GPD parameters exhibit some symmetry with respect to the middle of the domain? Could the lack of such symmetry be related to problems of statistical convergence for the estimators of these parameters?

- Fig. 7 caption: "from the long DNS (dashed black curves)": did you mean short DNS?

- p30, "displayed in Fig. 9a": the figure does not have multiple panels.

- section 5.2 and caption of Fig. 8: the notation R^2 for the coefficient of determination is a bit confusing given that R is the intensity, even if you never need to square it...

- Fig. 9: for long ASTs and small scale parameter $s$ the conditional severity PDFs have most of their mass outside of the sampled severities. Is this robust from a statistical point of view?

- Fig. 10: I am not fond of the unusual scale used for the two bottom rows, which at first glance gives the impression that the correlation coefficient depends essentially linearly on the AST.

- Fig. 11 caption: "Shaded regions show the areas between truncated upper and lower means." I am not sure I understand the justification for this.

General comments:

This manuscript studies the estimation of rare event probabilities in chaotic and high-dimensional systems, such as climate models. This has been a long-standing challenge: naive sampling methods can accurately estimate the probability of high-probability (or frequent) events, but they lead to large errors when estimating low-probability (i.e., rare) events.

Among existing methods, this manuscript focuses on rare event sampling (RES) with ensemble boosting. As the authors point out, two immediate questions arise: 1) What type of perturbations should one introduce to create the ensemble, yet avoid unrealistic rare events. 2) At what times should the perturbations be applied (Advanced Split Time, AST)? The authors mostly focus on the second question and propose/investigate three heuristics for choosing an optimal AST. Two of these methods are threshold based, whereas the third one is more systematic and selects AST by optimizing a prescribed functional. The authors propose two such functionals: Eq. (26), Expected Improvement (EI) and Eq. (27), Thresholded Entropy (TE). The proposed methods are investigated exhaustively on a two-layer quasi-geostrophic (2LQG) model. 

Overall, this manuscript addresses an important aspect of RES based on ensemble boosting. The proposed criteria are motivated well and seem sensible. Furthermore, the authors are admirably detailed in describing their methods and the related numerical results. However, I have a few comments/questions that will hopefully improve the manuscript.

Specific comments:

- In principle, the rare event probabilities can be estimated to an arbitrary accuracy if the DNS simulations are long enough. Of course, this would be computationally expensive for exceedingly rare events. Therefore, the main purpose of RES is to reduce the computational cost for a desired accuracy. With that in mind, I was surprised to find no data (table/figure) reporting the computational cost of different methods. The authors should report these in the manuscript since computational cost is a vital piece of information for this work. 

- Related to the previous comment: if I understand correctly, the equal-cost curves in Fig. 13 indicate that the ensemble boosted RES with the optimal AST is almost as good as a long DNS simulation. In fact, in this figure the RES results do not capture the most extreme events (severity>.67), which makes me think the long DNS is doing even better than RES. At any rate, if at the same computational cost, the DNS results provide us with the same information (or even more) as RES, why should one bother with RES and finding an optimal AST? 

- I understand that the type of perturbations (spatial structure) is not the focus of this paper, but it would be helpful to comment on its impact on AST. For example, do you expect the optimal AST to change if a different type of perturbation is used? Why or why not? The answer to this question is particularly relevant since the perturbations in section 4.1 arise from the most unstable mode around the baroclinically background state. This mode has no physical relevance to the chaotic trajectories being perturbed (since they are presumably far away from the background state). 

- In Section 3.3 the target variable is introduced as upper-level concentration averaged over a small spatial box. Does this target variable have some physical significance? If so, please explain it more clearly. If not, I'm wondering what the motivation was for choosing this particular ``intensity function of interest''.

- Although the authors do a good job reviewing the relevant literature, there are two prominent methods that are missing. One is a rare event estimation method based on LDT, see e.g., 

(a) G. Dematteis,T. Grafke, & E. Vanden-Eijnden,  Rogue waves and large deviations in deep sea, Proc. Natl. Acad. Sci. U.S.A. 115 (5) 855-860 (2018).

(b) Tong, Shanyin, Eric Vanden-Eijnden, and Georg Stadler. "Extreme event probability estimation using PDE-constrained optimization and large deviation theory, with application to tsunamis." Communications in Applied Mathematics and Computational Science 16.2 181-225 (2021).

Another related method is the variational approach introduced in

(a) Mohammad Farazmand, Themistoklis P. Sapsis ,A variational approach to probing extreme events in turbulent dynamical systems.Sci. Adv.3,e1701533 (2017)

(b) Blonigan, Patrick J., Mohammad Farazmand, and Themistoklis P. Sapsis. "Are extreme dissipation events predictable in turbulent fluid flows?." Physical Review Fluids 4.4 044606 (2019).

Both above methods seek to identify the onset of rare events in a systematic way and therefore closely relevant to the present manuscript. I don't mean for the authors to compare their RES against these methods, but some discussion seems appropriate.

Technical corrections:

- Page 14, Eqs (26-27): It is mentioned that the EI and TE are ``optimized''. Please clarify whether they are minimized or maximized.

- Page 15: ``We have conjectured that ...'' I believe the authors are conjecturing in this manuscript for the first time, but the word ``have'' makes it sound like it was conjectured in a previous study. If so, please add a reference.

- Page 16: ``The model setup aims to distill some challenges we have encountered with rare event algorithms across the hierarchy.'' It is unclear what you mean by ``across the hierarchy''.

- Eqs (28-33): A 2LQG model should have a stream function (and passive scalar density) for each layer. But this model seems to have only one stream function. Is there a subscript missing from the dependent variables? 

- On a more general note, I encourage the authors to edit the manuscript thoroughly. There are numerous sentences/paragraphs whose writing can be improved significantly and therefore improve the readability of the paper.