The manuscript under consideration is focused on an open problem, the prediction of large earthquakes. This paper presents an interesting study based on Machine Learning, ML, to predict intense earthquakes in terms of the released energy.

In  strictly sense, prediction is considered when the position of the epicenter, the time of occurrence and the magnitude are known with precision.

As it  is well known, the underlying dynamics of tectonic plates is due to complex processes inside the Earth such as convection. Those involved processes  give rise to the interaction between tectonic plates, being stick-slip the main mechanism for the  earthquakes occurrence. On the basis of these aspects, seismic dynamics is complex.

This paper shows a good prediction approximation produced by ML-based algorithms.

The work is very important and their results are very interesting, however, in my opinion, some important aspects have not been considered in this study:

a) The catalogue use for the study considers a period from 1900 to 2015. In the period from 1900 to 1920 very few earthquakes are observed while in recent years the density of events is much higher. Authors should explain these differences  time periods.

b) Figure 1a shows the magnitude of completeness (M6.4) which is only part of the Gutenberg-Richter law. The authors should show the Gutenberg-Richter law taking lower magnitudes (for example, M < 3) and thus determine the correct magnitude of completeness. To do this they can use the ZMAP platform where they can estimate the b-value and completeness Mc of the GR law.

c) It is important that the authors justify the reason for considering only earthquakes with large magnitudes and avoid the other ones with low magnitudes.

d) To calculate the IMF functions authors used the algorithm proposed by Huang et al. (1998). In my opinion, the authors should explain how the complex dynamics of seismic activity is involved to obtain IMF to obtain the smooth curve (a) in Figure 2.

e) When applying the algorithm proposed by Huang et al. it is not clear the role seismic activity with magnitudes M < Mc could have.

Their results do not consider time and epicentral position.

I consider that these suggestions can improve the context of the study that the authors have developed and then could be published

The manuscript considers the problem of seismic energy forecasting.

The method proposed builds on previously published research while improving the results.

Overall, the manuscript is well structured, although some improvements are needed to improve the readability and reproducibility of the results.

Section 1 (Introduction) and Section 2 (Background) cover, for the most part, the same topics: why are they divided?

Furthermore, while these sections may provide a partial overview of previous studies using ML for seismic energy forecasting, these may be presented in a more coherent form.

In Section 1, the concept of "ensemble" is introduced in a way such that it seems like a modern concept, despite being used in many fields for many decades.

At line 158, the authors mention the inclusion of expert opinions in the random forest (RF): the authors should specify if this is a general comment on the RF approach or for this study. If the latter is true, how has the expert opinion been included?

In section 3, do the authors consider the uncertainties due to poor station coverage, empirical relationships, etc.?

The authors should discuss how completely removing the events with magnitudes below the magnitude of completeness could improve the results: having information about lower magnitude earthquakes, even though it is partial, should still be better than not having any information at all. Also, this could improve the temporal distribution of the events, possibly allowing for monthly basis analysis.

In Table 3, how is the %(variance) being computed?

In Section 4, the authors should describe in further detail the inputs and outputs of the different models. As an example, the authors should explain in greater detail the role of the lag and how the inputs are computed, especially at the beginning of the timeseries. Furthermore, the authors should better explain what the inputs and outputs are during the testing phase.

In the MLP section, considering the few data points in the training dataset, doesn't a batch of 100 samples consist of the whole dataset (especially for the Western Himalaya case)?

In the Linear Regression section, the authors should add a reference to the M5 method. Furthermore, the ridge hyperparameter seems to be specific to the authors' implementation, rather than a general parameter of LR.

The parameters S, Z, and Y should be properly introduced at the beginning of the section for improved readability.

What does a_r represent in equation 5? Is it the feature vector?

Check the formatting of equation 6: there's a missing parenthesis.

In Section 5, the use of the training dataset to choose the best model should be carefully motivated. This practice can produce biased models as they are selected on the same data they are trained on, hence favouring overfitting models. A better solution would be to use a third independent dataset (i.e. the validation dataset).

Section 8 (Discussion) could be extended, considering how the proposed method could be further improved in light of its current limitations.