General comments:

The paper describes a significant step in the further development of the magnetospheric simulation model Vlasiator, by adding an ionospheric module, similar to the one used in the MHD-based model GUMICS. Magnetospheric simulation models are useful for describing and predicting the impact of solar wind on the earth and its atmosphere, and there is currently much interest in the further development of these models.

The paper is well written and well organized. I have only some minor comments.

We thank the reviewer for their helpful and insightful comments on the manuscript, which we have addressed below. Our answers are marked in bold text, like this line.

Specific comments:

Section 2:

Can you say something (very roughly) about the hardware requirements of Vlasiator, such as computer power and processing time? How feasible is it for the reader to try to run a simulation?

Ah, this is indeed an important stepping stone that still hampers larger-scale adoption of hybrid-Vlasov simulation methods in space physics: their large computational requirements.

As not to give anybody false hopes, a small sentence has been added, outlining that global simulations with Vlasiator are requiring some millions of core-hours computing time on supercomputers.

Figure 2:

The right-hand graph shows a refinement different from what you describe in the text (line 125-128). It may be that you mean this figure only as an example of what can be done. But it is confusing. It would be nice if the figure showed the exact refinement areas that you are actually using.

Figure 2 indeed shows a quite simplified mesh, with n=256 base grid points and its subsequent refinement stages, to illustrate the refinement process. We had indeed considered presenting the full mesh resolution of our global simulation ionosphere mesh in this figure, but found that it became very hard to discern any features at those resolutions. Hence we decided to retain this low-resolution example as Figure 2, and introduced the polar-cap closeup in Figure 3 in order to actually present the mesh geometry that is employed in our magnetosphere runs.

In order to reduce the confusion, the image caption of Figure 2 has been amended to explicitly reference Figure 3 for the full, production-scale resolution mesh.

Line 140: 

Where exactly is the upper endpoint of the field lines, to which you interpolate the magnetospheric quantities? Is it exactly at the point where the traced fieldline reaches r = rC?

Yes. This is a deliberate choice to be independent of the mesh geometry of the outer magnetospheric simulation. Magnetospheric properties at this location are then obtained by interpolation at said location.

Line 157:

What is the “downwards-facing hemisphere”? Do you mean: the hemisphere where the fieldlines point earthward?

Indeed, since the goal here is to formulate a model for the precipitating electron velocity distribution function shape, only the half-space of velocities that are moving towards the Earth are relevant for precipitation.

We replaced "downwards-facing hemisphere" by "the velocity half-space where particles move earthward", hoping that this formulation is clearer.

Line 215-219:

Please explain this reasoning more. If different sources of conductivity are dominating at different altitudes, why is it then better to add up the different contributions according to Eqs. (15)-(16) rather than just summing them linearly?

In our model, the free electron density n_e is assumed to be proportional to the square root of the ionoization rate (eq. 10). This value then linearly contributes to the calculation of the conductances (eqs. 11-13).

In the end, we wish to sum the contributions of different ionisation processes (starlight, EUV and particle precipitation), and it would make more sense to sum their individual ionisation production rates q. But alas, the semi-empirical models that we employ only provide for height-integrated conductivity contributions, which we assume to be similarly proportional to sqrt(q).

And even though we state that the ionisation effects of these processes are dominating in different atmospheric layers, we can unfortunately not do any better with the data that our modeling provides. Thus, we sum in quadrature.

While this might be a somewhat coarse approximation, this is the standard trick of the trade likewise employed in all other global simulation codes with height-integrated ionosphere models that we have encountered in our literature search (compare references Ridley et al., 2004, Janhunen et al. 2012).

Line 223-238:

You describe three approaches to place the conductance tensor in the geometric coordinate system. You describe problems associated with methods 2 and 3, and therefore choose to use method 1. However, you do not mention any possible downside of method 1. In particular: when B is in reality not vertical (especially at lower latitudes), perpendicular conductivities will be partly in the vertical direction (which part will presumably be ignored). Can you explain how you deal with this? Is maybe the equatorial part less important anyway, because the field lines do not reach as far out?

Since the ionosphere is treated in a height-integrated manner, the concept of a "horizontal B component" does not actually make sense in the model's context.  Essentially, the height integration is not happening in radial direction, but along the field line direction. But by *assuming* a radial B direction in this step, the details of this process are ignored and potential cos(theta) contributions to effective atmospheric height are neglected. This is assumed to be of no major detriment close to the auroral regions, as the magnetic field direction is still close to radial. The more equatorial regions, in which this might have an effect, are not contributing to the polar ionospheric dynamics in our model, as their magnetic field lines close onto themselves beneath the magnetospheric simulation inner boundary anyway (and are thus set to zero). In fact, the shielding latitude value Theta_Shield, discussed further down, is partially motivated by these considerations.

Line 254: 

Does your solution of putting a constraint only “alleviate” the instability? I would think that this solution (effectively: putting a boundary condition on an  underdetermined system) completely stabilizes the system.

Indeed, the choice of the word "alleviate" was unfortunate, it has been replaced by the word "fix".

Line 255:

“Theta_Shield in [0 ... 70]” do you mean that Theta_Shield can be chosen anywhere between 0° and 70°? Or that you mean to make the potential zero at all latitudes from 0° to 70°?

Yes, the shielding latitude can be chosen between 0° (effectively only pinning a small strip of mesh nodes on the equator to zero latitude) and 70°.

Higher shielding latitude levels are numerically more efficient (as a smaller fraction of the ionosphere mesh actually needs to be solved), but care must be taken not to clip away any of the actual interesting auroral regions.

Line 271:

“the first layer of cells adjacent to the ionospheric boundary” Do you mean: the cells adjacent to the inner boundary cells?

Indeed. It has been reformulated as "the magnetosphere's inner boundary".

Line 314:

“scatterplot.” -> “scatterplot, as a function of latitude.”

Great suggestion! The text has been adapted as such.

Line 328:

“Maxwellian velocity distribution” Velocity in which direction? (Presumably in the direction of the ExB drift?)

Since the simulation domain is fully 6-dimensional, the state of each spatial simulation cell needs to be initialized by a choice of 3D velocity distribution function. The Maxwellian distribution function in our case is chosen to be isotropically centered around the nominal solar wind speed, in -x direction in our coordinate frame. 

The text has been reformulated, to make it more explicitly clear that the solar wind velocity is taken as the centre of the Maxwellian VDF in every cell.

Line 374:

“…in a converse manner…” What you describe after this is not really a “converse” (opposite) manner. I suggest to write: “The mesh resolution N can be adjusted in order to resolve physical phenomena...”

Agreed, the term "converse" might be confusing here, and the sentence has been reformulated similar to your suggestion.

Line 378:

“The result in Figure 8 shows correctness”

“correctness” sounds very absolute. I suggest:

“The result in Figure 8 shows overall correctness”

The word "overall" was added, as you suggested.

Line 404-405:

“This is because the drift motions of ions and electrons, which should be oppositely directed, cannot currently be taken into account separately” How do you mean “should be”? Are they not oppositely directed in MHD simulations?

Also, I would think the problem is not so much that their directions are opposite, but mostly that their velocities can be different, which MHD simulation doesn’t account for. Isn’t that so?

In the typical single-fluid MHD approach, there is no distinction of different particle populations, their velocities or current contributions, as MHD only deals with overall mass- and momentum flows (and potentially higher moments).  As such, MHD can not model distinct dawn- and duskward drift motions of ions and electrons.

A kinetic method, like Vlasiator's, could in principle model these distinct drift motions. Our choice of a hybrid approach however, in which electron dynamics are reduced to that of a massless, charge-neutralizing fluid, only provides us with kinetic information of the ion species. For lack of any other option, we use the ion data as a proxy quantity for electron precipitation. We hence observe asymmetries in FAC and conductivity patterns that appear mirrored wrt. physical expectations, and interpret these as indicative of mishandled drift motions.

Line 433-434:

“Some preliminary evidence of new kinetic-physics features is present”

What do you mean? Did you present any evidence of this in the paper? If so, where? If you are referring to evidence which you have but are not showing, please either describe it a bit more, or otherwise simply mention that work is ongoing to develop this feature.



The sentence in question has been reformulated to clarify that analysis work is ongoing and will be topic of future publications.

Technical corrections:

Line 18:

“space weather” is written twice.

The duplication was removed.

Line 119:

Probably you mean Ri instead of RE here. (Although it won’t make much of a difference.)



Indeed! The formula has been corrected accordingly.

Line 130-134:

The word “trace” seems to be used with different meanings. In line 131 you

trace the node coordinates along the fieldlines, and in line 133 you trace the

fieldlines. Possibly in the sentence “ionosphere mesh node coordinates are

traced upwards along the fieldlines”, “traced” is not the right word.



The sentence in question has been reformulated to clarify that the node

coordinates are the starting points for stepping along the field lines.

Line 149:

“electron distribution” -> “electron velocity distribution”



This clarification has been added,

Equation (3):

The k of the first kB should probably also be in italic.



Indeed. It has been reformatted such.

Line 255:

“[0o . . .70circ]”

Something seemed to go wrong with the degree sign after the “70”.



This formatting mistake has been fixed.

Line 289:

remove “(up to a proportionality constant)”

When you wrote “the same shape”, you already implied this; besides, Φ and j do not even have the same unit, so there must necessarily be a factor between them.



Agreed, the clause in parentheses was redundant and has been removed.

Equation (24):

The J of J0 should probably be in lowercase.



Indeed, that is a reasonable choice to keep the terminology consistent. It has been corrected.

Line 314:

“Φ a” -> “Φ in a”



Oh, yes indeed. Fixed.

Line 332:

“rC = 35.7km”

Presumably, you mean 35.7x103 km?



Indeed, the value has now been scaled correctly.

Line 357:

”increase both” -> “increase in both”



Corrected.

Line 371:

Presumably, you mean with the term “multipole” the harmonics with l and m > 2

as shown in figure 6. But you didn’t use that term in Section 3.1, so it’s a bit confusing to refer to it now with that term.



The term "multipole" has been replaced by "spherical harmonic" in all instances that it occurred in the document. This should now be consistent throughout the text.

Line 385:

“without” -> “with” (you’re saying already “=0”)



That makes sense! The text was corrected accordingly.

Line 385:

Capitalize “H” in “Hall”



Corrected.

Line 407-408:

You are using the word “model” twice with different meanings in this sentence.

I suggest:

“Work is ongoing to develop a more sophisticated precipitation model that takes the plethora of kinetic simulation data from Vlasiator more effectively into account and to provide for a more realistic precipitation distribution function (e.g., Zhang et al., 2015).”



Excellent, the formulation was adopted as you suggested.

Line 409:

“Further research will be required to overcome this limitation.”

This sentence seems to say the same as the previous sentence, and can probably

be removed.



The sentence has been removed.

Line 420:

The sentence starting with “In the meantime…” talks about a different subject, so it would be good to start a new paragraph.

Actually, it seems this refers to the same subject as the paragraph of lines 399-411 (kinetic proton information), so maybe it belongs there.



The paragraphs have been slightly reordered, with the sentence in question now directly following the discussion of precipitation model developments.

General Comments:

This paper describes and performs basic verification for the first Vlasiator magnetosphere-ionosphere coupled simulation. Drawing from GUMICS heritage, the ionosphere potential solver is built from scratch. The model is described in detail, and two verification tests show convergence and agreement with known solutions. The complete magnetosphere-ionosphere simulation is then run for 500s, showing development of expected current systems.

This paper is thorough, well-written, and complete. It is a critical documentation of new Vlasiator capabilities at early steps, and will form an important reference for science investigations to come.



We thank the reviewer for their helpful and insightful comments on the manuscript, which we have addressed below. Our answers are typeset in bold, like this line here.

 

Specific Comments:

[Lines 77-89]:

More citations needed? Also, this is entirely missing any mention of the work of the Michigan group, who are currently doing coupled MHD and PIC simulations as well as advanced ionosphere conductance modeling.  Additional literature review would be useful both here and above to establish the state of the whole field, not just the GUMICS heritage.



Thank you for the suggestion. While we did include a reference to the Michigan group's SWMF framework, we agree that this is insufficient, especially in light of their extensive work on coupled MHD / PIC simulations, which can also interact with ionospheric model implementations (although they are never directly coupling their PIC domains with the ionospheric model). We have added a sentence highlighting these modelling efforts and putting them into context with the presented new model.

Furthermore, a sentence has been added that further highlights Michigan's SWMF as well as APL's MAGE as the state of the art in terms of magnetosphere<->ionosphere coupling.

[Line 115]:

Would you ever choose 2^24 mesh points? What is the purpose of this statement?



The listed limit of 2^24 mesh points had been determined as the maximum numerical stability threshold of the presented spherical Fibonacci discretization method. At mesh resolutions beyond it, the listed properties break down, as presented in the Keinert et al. (2015) reference.

An actual simulation mesh of this size would be completely impractical due to its immense memory and computational requirements.

The text has been updated to clarify that this number defines a theoretical upper limit.

[Line 255]:

What are the impacts of choosing different latitudes in this range? Can you offer guidance to make this choice?



A more poleward shielding latitude increases numerical efficiency of the solver, as a smaller domain needs to be solved, but care must be taken not to cut off current closure in any physically relevant regions.

Due to the Magnetospheric simulation's inner boundary being situated at ~4.7 R_E, the field lines at low latitudes are not connecting to the Vlasov simulation at all, but are closed beneath it. Hence, their field-aligned current is set to zero, and the ionospheric mesh below 60 degrees does not contribute to auroral dynamics. Hence, it makes sense to chose the shielding latitude a bit equatorward of this boundary.

[Line 264]:

This binary choice of corotation is interesting. Also, why is this L=5 different from the coupling radius of 5.6 Re?



The plasmapause, marking the spatial limit of corotating plasmasphere populations around Earth, indeed presents a rather sudden discontinuity in plasma bulk velocity. While detailed models of its electric model exist (such as the one provided by the Maus et al. (2017) paper given in the text), the simplified choice of a hard limit at L=5 was made to provide at least some effect of the corotation into the magnetospheric simulation box.

The coupling radius (forming a spherical shell at 5.6 Re) is conceptually disconnected from the toroidal shape of the plasmasphere region. The choice of coupling radius is predominantly motivated by numerical stability concerns of the bidirectionally coupled ionosphere<->magnetosphere system.

```

[Line 335]:

More mention of Figure 10? What is the purpose of including it? Right now, it is only briefly described in Line 335.



The purpose of Figure 10 is to illustrate the complete set of input quantities that are downmapped from the magnetospheric domain and used by the ionospheric

model to solve the potential distribution. You are correct that these quantities' phenomenologies are not further discussed in this manuscript, which focuses mostly on the methodical aspects of the solver.

[Line 344]:

What is the clear improvement? Can you provide a figure for comparison, instead of just a reference? This should be a major point of this paper.



We agree that the previous formulation was not clearly presenting the qualitative improvement of Vlasiator's reproduction of ionospheric current

systems. A more specific reference to Figure 3 of Horaites et al. (2023) has been added.

While we agree that a comparison figure might help the reader to directly appraise the change with the implementation of the new model, we have decided

against it in the presented manuscript, as its focus lies fully on the properties of the new method.

Note also the paragraph on the current systems in the discussion section of the manuscript (lines 398-408), which further elaborate on their dynamic evolution.

[Section 3.3.4]:

This is a very short section. Anything else to say here?



This section has been deliberately held short, as the presented manuscript's focus lies on the model description and validation. Hence, only the

initialization phase of the global simulation presented here is discussed, until an initial steady state has formed.

The dynamic evolution of magnetospheric phenomena, and their interaction with the ionosphere will be the content of multiple upcoming presentations,

focusing specifically on their physical contents. We have added a short paragraph outlining this fact.

[Line 433]:

What evidence? This is not mentioned earlier, as far as I can tell. Please be specific, and add this somewhere, or remove the mention.



The sentence in question has been reformulated to clarify that analysis work is ongoing and will be topic of future publications.Technical Corrections:

[Line 18]: The words ‘space weather’ are repeated.



The duplication was removed.

[Line 130-135, elsewhere]: Should ‘fieldlines’ be ‘field lines’?



All occurrences have now been replaced by "field lines", to have a consistent terminology in the manuscript.

[Line 148]:

Consider changing ‘electron density’ to ‘electron number density’ for clarity.



Excellent suggestion, this has been added.

[Line 175]:

As a general rule, parenthetical expressions can be incorporated into actual sentences. In this line, the parentheticals make the statement difficult to

read. Consider splitting this into multiple sentences that incorporate the parenthetical information into the main points of the sentences.



The parenthetical clause has been merged into the sentence, and the resulting long sentence split in two.

[Line 180]: Please change the colon to a period, since this is a paragraph break.



Understood, it has been changed.

[Line 234]: Remove the periods after 2 and 3.



Ok, they have been reformatted without periods.

[Line 255]: The degree symbol near 70 appears to be improperly coded.



This formatting mistake has been fixed.

[Line 278]: Replace the colon with a period.



Ok, it has been replaced.

[Line 295]: The parentheses are unnecessary.



They have been removed.

[Line 332]: There seems to be an error in the coupling radius in km.



Indeed, the value has now been scaled correctly.

[Line 357]:

Missing the word ‘in’. Also, consider calling those features ‘local peaks’ or ‘local maxima’ instead of ‘slight increases,’ since the whole sector is not slightly increased. Only the local times highlighted experience the local peaks.



The missing "in" has been added, and the text reformatted per your suggestion.

[Line 373]: Remove capital ‘I’ in ‘Ionospheric’ – this is not a proper noun.



The spelling has been corrected.

[Line 383]: Remove the second comma.



The comma was removed.

[Line 385]: Capitalize Hall.



Corrected.