We thank the referee for their thoughts and inputs. Some revisions were done on that basis, attempting to increase the quality of the paper:

1. The initialisation of minor neutral constituents is clarified and the NRLMSIS2.1 model is used to include NO. The sensitivity to uncertain initial NO density is tested. Minor neutral constituents not represented in NRLMSIS2.1 are initialised with 0 density, with a leading time of half an hour to let them approach steady-state densities. This may not be long enough for all species to reach steady-state, but it also not a requirement. The aim is to construct a reasonable close model of the ionosphere, in regard of uncertain initial conditions. In the auroral region, where precipitation varies spatially and changes rapidly with time, the production of minor species will have a corresponding spatial and temporal variation. Therefore, we cannot expect all species to be at steady-state.

We tested the robustness of this approach towards uncertain initial conditions, and found that the differences are very small.

2. Reaching equilibrium (or steady-state) is not a requirement. While we have uncertain initial conditions and only initialise densities that we have models for, the leading time is intended to mitigate gross mismatches that would otherwise lead to large corrections at the beginning of the dataset, and letting the uninitialised species reach non-zero densities. The largest corrections should occur at the start. Furthermore, both NO and N(4S) are theorised to never reach equilibrium in the limited time between substorms. Auroral precipitation may increase NO and N(4S) concentrations considerably (Bailey et al, 2002). Lastly, due to the rapid variations in precipitation, many other species will also not be in equilibrium over the entire evaluation period. For all these reasons, we feel that our approach of giving the model a limited time to approach steady-state conditions is justified, without the need of actually being at steady-state. This ensures that the all species densities are in a reasonable range.

On the time-scales of auroral precipitation, vertical diffusion is small at E-region heights, see e.g. Turunen et al. (2009).

Bernhardt et al. (2000) showed that an artificial airglow is subject to neutral diffusion in all directions, and drifts in response to neutral winds. Ions and electrons are subject to electrodynamic forces present in active aurora (e.g. Gustavsson et al. 2001, Krcelic et al. 2014). Furthermore, there are steep horizontal gradients in auroral precipitation. Therefore we see it more useful to include diffusion, drift and convection in 3 dimensions. This will be incorporated in future work, depending on EISCAT3D data availability.

3. The reaction rates are updated. Future work will include more reactions, in particular decay of excited states.

References:

Krcelic, P., Fear, R. C., Whiter, D., Lanchester, B., & Brindley, N. (2024). Variability in the electrodynamics of the small scale auroral arc. Journal of Geophysical Research: Space Physics, 129, e2024JA032623. https://doi.org/10.1029/2024JA032623

Gustavsson, B., et al. (2001), First tomographic estimate of volume distribution of HF-pump enhanced airglow emission, J. Geophys. Res., 106(A12), 29105–29123, doi:10.1029/2000JA900167.

Bernhardt, P. A., M. Wong, J. D. Huba, B. G. Fejer, L. S. Wagner, J. A. Goldstein, C. A. Selcher, V. L. Frolov, and E. N. Sergeev (2000), Optical remote sensing of the thermosphere with HF pumped artificial airglow, J. Geophys. Res., 105(A5), 10657–10671, doi:10.1029/1999JA000366.

Bailey, S. M., C. A. Barth, and S. C. Solomon, A model of nitric oxide in the lower thermosphere, J. Geophys. Res., 107(A8), doi:10.1029/2001JA000258, 2002.

E. Turunen et al., “Impact of different energies of precipitating particles on NOx generation in the middle and upper atmosphere during geomagnetic storms,” Journal of Atmospheric and Solar-Terrestrial Physics, vol. 71, no. 10, pp. 1176–1189, July 2009, doi: 10.1016/j.jastp.2008.07.005.

We thank the referee for their insightful and helpful comments and suggestions.

1. The ElSpec methods are not based off of standard inverse-problem regularisation techniques

like Tihkonov regularisation or Maximimum Entropy, but instead uses Akaike Information

Criteria for model-selection to handle the ill-posedness of the estimation problem and produce

simple (in the Occam's razor sense) electron spectra. This term is not a regularisation term, but

a penalty that prevents overfitting. This corrected penalty term accurately accounts for small

samples sizes, . It is this sample-size correction that distinguishes the corrected AIC, cAIC, from

the standard AIC. For the detailed description see the first publication on ELSPEC (Virtanen et al., 2018),

or Burnham et al. (2002): Model Selection and Multimodel Inference (2nd ed.).

2. The entire objective with this figure is to show that the different initial conditions are producing

almost exactly identical densities, such that there are no discernible differences in Figure 4. We will clarify this point in the final version.

3. Percent differences will blow up when the reference spectrum goes to 0, i.e. at the high- and low-energy tail. We added time slices of the absolute fluxes and the difference, such that the relative difference can be better seen. 

4. The focus of this paper is developing new methods for improving the inversion of

electron density to electron spectra during aurora, incorporating the E-region chemistry.

The improvement might be considered minor, but even if it is a small

improvement it is a significant improvement, leading to a systematic correction

of the estimated primary electron spectra. Every simplifying approximation made

has to be scientifically justified, and then be cautiously kept in mind throughout its

application with careful checks that the conditions of the approximation is fulfilled. We

prefer to avoid this approximation and the ensuing worries that follow. As such we argue

that this endeavour is worth to publish, and to our knowledge, has not been done before.

Other authors have taken variable effective recombination rates into account (see e.g.

McLennan et al.), but not by using fully dynamic model of E-region chemistry. The

mentioned GPI model (Kaeppler et al.) uses a 5 species model, seems, to our understanding,

to be more tailored to D-region studies, since it models the densities of positive and negative

 "light" and "heavy, cluster" ions. Such a model would be ill-suited to resolve the dynamics of

NO+ and O2+, which determine the recombination rate in the E-region.

 E-region chemistry on its own, on the other hand, which is the context of Figure 10, we

 agree, has been done before. Therefore Figure 10 is not the focus of this paper, albeit

 an interesting finding.

5. See Point 4.

References:

J. McLennan, A. Jaynes, R. Troyer, S. Kaeppler and M. Shumko, "Validating an Energy Flux Inversion Method with Satellite Data (Invited)," 2025 URSI Asia-Pacific Radio Science Meeting (AP-RASC), Sydney, Australia, 2025, pp. 1-4, doi: 10.46620/URSIAPRASC25/YWNE9813.

S. R. Kaeppler, R. Marshall, E. R. Sanchez, D. H. Juarez Madera, R. Troyer, and A. N. Jaynes, “pyGPI5: A python D- and E-region chemistry and ionization model,” Front. Astron. Space Sci., vol. 9, Dec. 2022, doi: 10.3389/fspas.2022.1028042.

I. I. Virtanen, B. Gustavsson, A. Aikio, A. Kero, K. Asamura, and Y. Ogawa, “Electron Energy Spectrum and Auroral Power Estimation From Incoherent Scatter Radar Measurements,” Journal of Geophysical Research: Space Physics, vol. 123, no. 8, pp. 6865–6887, 2018, doi: 10.1029/2018JA025636.

Burnham, K. P., & Anderson, D. R. (2002). Model Selection and Multimodel Inference (2nd ed.). New York: Springer.