The manuscript “Reconstruction of Mercury’s internal magnetic field beyond the octupole” by Toepfer et al. presents and compares different existing inversion techniques to reconstruct the internal magnetic field data of Mercury. The authors use a simulated magnetospheric model to get the synthetic magnetic field data for the inversions, to mimic the BepiColombo (MPO) mission data acquisition in a non current free environment. In this way, the authors are capable of evaluating the different methods in retrieving the known coefficients given a priori to the simulations. This study is of great importance specially that BepiColombo is on its way to Mercury. Besides some concerns detailed below, the manuscript is well organized and well written. 

Main comments:

My main comment concerns the application to simulated Mercury magnetic field data. The authors make use of known Gauss coefficients from literature (Anderson et al. 2012, Thébault et al. 2018 and Wardinski et al. 2019), where the different models represent the Hermean internal magnetic field with Gauss coefficients of m=0. Those models are evidently constrained by the available data over a single hemisphere. 

As the authors are doing a theoretical exercise using a synthetic model, I wonder if the authors shouldn’t test the different inversion techniques for a more complex internal magnetic field, where m!=0 coefficients are also considered. For example, I have in mind the Jupiter’s internal magnetic field case where the magnetic field from a hemisphere is axisymmetric but not in the other. In the case Mother Nature surprises us with a more exotic Hermean internal magnetic field, are the different inversion methods described in this manuscript capable of representing the non-axial Gauss coefficients? I suggest to run a couple of simulations using non-axisymmetric coefficients to test the limits of the different inversion techniques combined with the Gauss-Mie parametrization.

My second main comment is related to the effectiveness of the different techniques to retrieve the Gauss Coefficients. This issue is actually highlighted when the authors used two different values for coefficient g^0_5. Depending on the strength of this coefficient the authors show that the inversion method perform differently. It should be mentioned that the used models (Anderson et al. 2012, Thébault et al. 2018 and Wardinski et al. 2019) give Gauss coefficients with a strong covariance, because there are no data available in the Southern Hemisphere. This means that several sets of coefficients are still valid to represent the internal field. Given that, I wonder what is the impact in changing slightly the many coefficient values (as done for g^0_5). How well are the different inversion techniques retrieving the given coefficient, and how much that value should vary to notice that change?  I suggest the authors to try other g_1^0, g_2^0, g_3^0, … values not too different from those already used, and that can potentially be also a solution for Mercury’s magnetic field, to check if there are limitations in discerning one from the other.

Moderate comments:

Introduction: I find that a description of the state of the art of the existing internal magnetic field models using MESSENGER data is missing. The lack of an introduction on this topic is enhanced later in section 4.1. I would also specify the limitations of those models, for example, the coefficients that are covarying.

Lines 96 - 100 Please, add a sentence to describe how the thin shell approximation is affecting the results.

Section 4.2 After running the simulations but before selecting the data for inversions, is there a noise added to each data point? What is the error level? Is there more than a value considered? How it affects the inversions? There is a brief mention on this topic at lines 368-371 but it is not satisfying.

Lines 309 - 312 It might help display a figure with the spherical grid points used for the inversions. Please, indicate the grid resolution used. Also, how much time shall we wait for BepiColombo to acquire enough data for your inversions?

Conclusion: An important aspect on internal field modeling is the time variation (or secular variation). This should also be mentioned here, even if this is not the scope of this manuscript. 

Minor comments:

Lines 14-17 there is a lack of citations.

Line 26 Wardinski et al 2019 also estimate the size of the core.

Lines 34-36 There is also a disadvantage, the authors are modeling the external sources that are given by the simulations only, which could not be a full representation of the real magnetosphere currents.

Lines 87 - 89 Please, define O in the text from equations 8 and 9.

Table 1 and 2: You have described in section 3.1 the Least Square Fit method, I would also add a column for that method in the given tables for comparison purposes.

Typos:

line 176 Imagenary -> Imaginary

The manuscript "Reconstruction of Mercury’s internal magnetic field beyond the octupole" by S. Toepfer and co-authors provides an interesting comparison between different inversion methods to reconstruct the internal magnetic field of Mercury. They used a hybrid plasma code to simulate Mercury’s magnetosphere, then they applied different methods to evaluate the "a priori" known coefficients. I think that the manuscript is a valuable contribution, especially for the BepiColombo community but not only. Indeed, these findings could be widely expanded to different environments, as well as, a similar methodological framework is a valuable support to any planetary mission. Furthermore, the manuscript is clear, well written, well organized, and well posed in terms of the existing literature. I would only suggest some improvements for the benefit of the reader as well as to assess their results.

1. I would recommend the authors to add errors on their coefficients' estimations as reported in Tables 1 and 2. Since both tables contains the main results of the paper, adding errors could improve the clarity of the results.
2. I would suggest the authors to comment on confidence intervals changes in their estimations under different solar wind conditions. Are they related to some specific solar wind parameters, apart the interplanetary magnetic field? I was wondering on plasma parameters like the Mach number(s) or the plasma beta. It would be nice to estimate coefficients under two/three different solar wind conditions or to add a few lines on this aspect.
3. Another possible interesting aspect to be mentioned could be the role of considering different harmonics degrees in terms of both l and m. Could the authors comment on the expected changes as a function of m and l?
4. I would suggest to add a few details on some specifics on the inversion models (noise, grid size, resolutions, ...) for clarity.
5. In light of the application of inversion methods to BepiColombo data, I would ask the authors to comment on the following aspect. The authors will use the model on magnetic field time series when MPO will explore different regions on the Hermean environment. I was wondering how the fast temporal variability of the Hermean environment as well as that of the different regions could affect the inversion methods used. Could this be considered as a "noise" for the method? Should be useful to firstly apply some filtering procedures on magnetic field data to remove the short-term variability and then apply the inversion method to the large-scale variability of MAG measurements? Could be some mixing between temporal and spatial scales that could affect the model performances?