Response to Ekaterina Svechnikova (CC1)

We sincerely thank Dr. Svechnikova for the careful assessment of the manuscript and the constructive suggestions.

Comment:

“The main limitation is that the degree of universality remains insufficiently demonstrated. The results are based on four high-altitude station configurations, one simulation framework, idealized vertically uniform electric-field layers, and a fixed post-field propagation distance. These assumptions are reasonable for a controlled scaling study, but they are much simpler than real thunderstorm environments, where electric fields are nonuniform, time-dependent, and often organized in complex multi-cell structures.”

Response:

The controlled and idealized configuration used in the present study was chosen to systematically examine the effects of atmospheric density and electric-field strength on RREA development and subsequent particle transport under comparable simulation conditions. The uniform vertical electric-field layers are not intended to reproduce the complete, event-specific electric-field structure of real thunderstorms. Rather, this configuration provides a controlled modeling framework because the spatial structure and vertical extent of strong electric fields inside thunderclouds are not sufficiently constrained by direct measurements.

The applicability of the results obtained within this framework should therefore not be extended to arbitrary thunderstorm configurations without additional validation. This limitation is now stated more explicitly in the revised manuscript, and the conclusions are formulated within the scope of the present CORSIKA simulation framework.

Comment:

“Several additions would considerably strengthen the result: comparison with other Monte Carlo tools, more explicit uncertainty estimates for the fitted parameters, sensitivity tests to field geometry and fitting range, analysis of different post-field gap lengths, and validation against observed TGE or gamma-ray glow spectra.”

Response:

The uncertainties of the fitted parameters are now explicitly reported, including the station-specific K coefficients and the parameters of the joint density-dependent fit, K=(6.39±0.26)×10³ kV and a=−0.923±0.083. In addition, the post-field transport analysis has been extended from the original 100 m separation to 25, 50, 100, and 200 m, allowing the distance dependence of the electron and gamma-ray spectra and particle counts to be examined directly.

A systematic comparison with independent Monte Carlo frameworks and sensitivity tests using nonuniform electric-field geometries would require additional simulations beyond the present CORSIKA configuration. Similarly, a quantitative comparison with observed TGE or gamma-ray glow spectra would require event-specific information on the thunderstorm electric-field structure, field termination height, atmospheric conditions, and detector response. These analyses are therefore not included in the present study.

Comment:

“Such checks would make the proposed scaling more robust, better calibrated, and more clearly positioned relative to existing RREA parameterizations.”

Response:

The revised manuscript places the proposed empirical scaling more explicitly in the context of previous RREA parameterizations and limits its interpretation to the atmospheric conditions and simulation configurations investigated here.

Comment:

“A minor stylistic point is that the manuscript is written by a single author but frequently uses collective formulations such as ‘we show’ and ‘we develop.’ Replacing these phrases with passive constructions or formulations such as ‘this work shows’ would improve stylistic clarity.”

Response:

The manuscript has been revised throughout to replace collective first-person formulations with passive constructions or neutral formulations, improving stylistic consistency.

 We sincerely thank the referee for the careful reading of the manuscript, the positive assessment, and the constructive suggestions.

Reviewer comment:

“Introduction, around the lines 26–27, recently discovered flickering gamma ray flashes (FGFs) should be added into the list of phenomena for completeness (e.g., Østgaard et al., 2024; Marisaldi et al., 2024).”

Response:

Flickering gamma-ray flashes (FGFs) have been added to the list of high-energy atmospheric phenomena in the Introduction, together with the suggested references to Østgaard et al. (2024) and Marisaldi et al. (2024).

Reviewer comment:

“It would be useful to briefly mention/discuss the fact that there are theories which consider more complex RREA, the so-called relativistic feedback mechanism (e.g., Dwyer, 2003).”

Response:

A brief discussion of the relativistic feedback mechanism has been added to the Introduction, including its role in initiating successive avalanches through backward-propagating positrons and high-energy photons, with reference to Dwyer (2003). The revised text also distinguishes this mechanism from the conventional RREA development considered in the present study.

Reviewer comment:

“Simulations are performed for four selected high-altitude observatories/stations. I recommend providing 1 or 2 references related to the observations of thunderstorm ground enhancement at each station.”

Response:

References to thunderstorm-related particle-enhancement observations have been added for all four high-altitude sites considered in the simulations: Aragats (Chilingarian et al., 2022, 2024), Nor Amberd (Chilingarian et al., 2025), Lomnický Štít (Chum et al., 2020; Kisvárdai et al., 2025), and LHAASO (LHAASO Collaboration, 2023).

Reviewer comment:

“Around lines 86–88, this simplified configuration should be justified by a reference.

Is the electric field (in Section 3) really constant/uniform within a vertical range of 2000 m in the simulations?

In addition, I believe that the range of 2000 m should be shortly discussed with respect to the avalanche length.”

Response:

The simulation configuration has been clarified in the revised manuscript. The electric field is modeled as a vertically oriented, uniform field throughout the 2000 m acceleration layer. The same idealized configuration was used in the previous CORSIKA study of Chilingarian et al. (2026), which is now cited in support of the adopted setup.

The role of the 2000 m field-layer thickness has also been clarified in relation to the simulated avalanche lengths. As shown in Sect. 3, the 2000 m layer is substantially longer than the effective avalanche lengths obtained for the electric-field conditions considered, allowing avalanche multiplication to develop within the acceleration region before the particles leave the field.

Reviewer comment:

“Equation on line 135. Introduce/specify H_st in the equation.”

Response:

H_st​ is now explicitly defined when it is first introduced.

Reviewer comment:

“Equation on line 178, specify i (which values of electric field and how many values were used?).”

Response:

The index i is now explicitly defined as labeling the four electric-field values used for each station. The values are 180, 190, 200, and 210 kV/m for LHAASO; 200, 210, 220, and 230 kV/m for Aragats; 210, 220, 230, and 240 kV/m for Lomnický Štít; and 230, 240, 250, and 260 kV/m for Nor Amberd. These values are now specified in the revised manuscript.

Reviewer comment:

“Sections 3 and 4, I suggest numbering all the equations that are written on separate lines (not only those, which are referenced later in the text).”

Response:

All equations displayed on separate lines throughout the manuscript have been numbered consistently in the revised version.

Reviewer comment:

“I think it would be useful to briefly discuss/show the change of energy spectra and counts for electrons and photons as a function of distance between the exit and detector (not only for 100 m).”

Response:

The post-field transport analysis has been extended to field-to-detector distances of 25, 50, 100, and 200 m. The revised manuscript now presents the evolution of both electron and gamma-ray counts and their energy spectra as a function of propagation distance beyond the lower boundary of the electric-field region. These additional simulations allow the different attenuation behavior of the electron and gamma-ray components to be systematically studied and extend the original analysis beyond the previously considered 100 m separation.   Minor and formal comments

Reviewer comment:

“line 41, … development of the relativistic runaway electron avalanche (RREA,... → …development of the relativistic runaway electron avalanche (RREA) theory…”

Response:

The wording has been corrected as suggested to “the development of the relativistic runaway electron avalanche (RREA) theory.”

Reviewer comment:

“Equation (1) It would useful to mention already here that E_th(h) depends on the air density.”

Response:

The density dependence of E_th​(h) is now stated explicitly when Eq. (1) is introduced, with

E_th​(h)=E₀​ρ(h)/ρ₀​.

We sincerely thank the referee for the careful reading of the manuscript and for raising important questions concerning the methodology and interpretation of the study. The comments helped us identify several points that required more precise explanation, and these have been clarified in the revised manuscript while preserving the scope and the main simulation results of the study.

Reviewer comment:

“The methodology suffers from a significant flaw. As I understand the manuscript, the parameterization used to characterize the field excess (as termed by the author) is not scaled consistently. Dwyer's formula (eq. 1) is a local relation that applies to reduced electric fields. In the present study, however, it appears that only the threshold field, Eth, is scaled with altitude, while the applied electric field is not.”

Response:

The electric-field ranges were selected separately for each observational site to obtain sustained supercritical RREA multiplication under the corresponding atmospheric conditions. This selection was guided by the previous CORSIKA study (Chilingarian et al., 2026), in which the RREA initiation threshold was investigated as a function of altitude and atmospheric density. In the present study, the runaway threshold is defined as

E_th(h)=E₀ · n(h),

where n(h)=ρ(h)/ρ₀ and E₀=284 kV/m. Thus, the density entering the empirical parameterization is a station-specific reference value․

The ranges used were 180–210 kV/m for LHAASO, 200–230 kV/m for Aragats, 210–240 kV/m for Lomnický Štít, and 230–260 kV/m for Nor Amberd. The revised manuscript now states this procedure explicitly.

Reviewer comment:

“The manuscript questions the invariance of the coefficient K, yet no physical justification is provided for why this invariance should fail. As written, the paper effectively claims a breakdown of similarity laws for RREAs. Such a claim requires compelling theoretical or numerical evidence.”

Response:

Applying the classical parameterization separately to the four stations gives K_station values ranging from 8.6×10³ to 11.2×10³ kV, with the station-specific fit uncertainties now reported in the revised manuscript. A simultaneous fit of Eq. (1) to all 16 simulation points gives K_fit≈9.58×10³ kV, which is used as reference for comparison with the generalized empirical relation.

The generalized relation was fitted jointly to the combined dataset, yielding K ≈ (6.39±0.26)×10³ kV and a = −0.923±0.083.

The factor n^a represents an additional density dependence beyond that already included through E_th = E₀ ·n.

Within the present CORSIKA dataset, this term reduces the inter-station scatter and increases the coefficient of determination from R²=0.90 to R²=0.99. Equation (10) is therefore presented as an empirical parameterization of the effective avalanche lengths obtained for the finite 2000 m field regions considered in this study․

Reviewer comment:

“Finally, a substantial body of relevant literature is omitted from the bibliography and discussion. The manuscript does not adequately place its results in the context of previous work.”

Response:

The literature review has been expanded to include previous studies of RREA microphysics, relativistic feedback, and Monte Carlo modeling. The discussion of the avalanche-length parameterization now includes Coleman and Dwyer (2006), including their density scaling and the parameterizations reported for different electric-field ranges. The revised Introduction also clarifies the applicability of the avalanche-length parameterization with respect to electric-field range and atmospheric conditions.

Specific and technical comments

Reviewer comment:

“l. 13. ‘Coefficient K’ is not defined at this point, nor is it defined later in the abstract.”

Response:

The Abstract has been revised so that the coefficient K is no longer introduced there before being defined. K is defined together with Eq. (1) in the Introduction.

Reviewer comment:

“l. 24. Replace ‘climatology’ with ‘meteorology’.”

Response:

This has been corrected.

Reviewer comment:

“l. 31. Not all of the listed phenomena are bursts; for example, gamma-ray glows are not.”

Response:

The wording has been revised so that the listed high-energy atmospheric phenomena are not collectively characterized as bursts.

Reviewer comment:

“l. 36. In addition to empirical models, many first-principles models (both self-consistent and non-self-consistent) have been developed and published.”

Response:

The corresponding statement has been revised, and the discussion has been expanded to acknowledge both empirical and first-principles approaches.

Reviewer comment:

“Equation (1) (and similarly throughout the manuscript): insert a space between the equation and the equation number.”

Response:

 The equation formatting has been corrected consistently throughout the manuscript.

Reviewer comment:

“l. 53. The most commonly used value of K in the literature is not simply 7200 kV, but probably that reported by Coleman and Dwyer (GRL, 33, L11810, 2006), who reported two different values depending on the electric-field range considered.”

Response:

 Coleman and Dwyer (2006) has been included and the discussion of the literature values of K has been revised accordingly.

Reviewer comment:

“l. 61. The complexity and variability of atmospheric electric fields are not, in themselves, an obstacle [...] What currently limits more accurate modeling and interpretation of observations is the lack of reliable measurements of the electric-field structure within thunderstorms.”

Response:

The statement has been revised to identify uncertainty in the thunderstorm electric-field structure as the principal limitation.

Reviewer comment:

“l. 63. Monte Carlo simulations have been extensively used to study RREAs and related phenomena in a large body of published work, none of which is cited here.”

Response:

The literature review has been expanded to include representative previous Monte Carlo studies of RREAs and related high-energy atmospheric phenomena.

Reviewer comment:

“l. 135. H_st is not defined.”

Response:

 H_st is now explicitly defined when first introduced.

Reviewer comment:

“l. 160. The term ‘extensively’ should be supported by appropriate references.”

Response:

Appropriate references have been added.

Reviewer comment:

“l. 163. The statement that ‘the universality of the coefficient K’ remains unresolved is misleading.”

Response:

The original wording has been revised. The manuscript now presents the analysis as an assessment of how consistently the classical empirical avalanche-length parameterization describes the effective avalanche lengths obtained from the CORSIKA simulations across the atmospheric conditions considered in this study.

Reviewer comment:

“l. 199. Please compare these results with previous studies.”

Response:

The revised manuscript now includes an explicit comparison with previous RREA calculations using the conventional density-scaled variables E_z/n and n·λ_eff. Representative results are compared with the ranges reported in Coleman and Dwyer (2006), Dwyer et al. (2012), and Skeltved et al. (2014).

Reviewer comment:

“l. 200. This discrepancy arises because the formula is not scaled consistently.”

Response:

This point is addressed in the response to the first Major Comment and has been clarified accordingly in the revised manuscript.

Reviewer comment:

“l. 284. E_e^max  should be defined before being introduced.”

Response:

 E^max, including E_e^max  and E_γ^max , is now explicitly defined before use as the characteristic maximum energy determined from the spectral distributions. In the revised manuscript, it is taken as the highest energy for which at least three consecutive energy bins contain more than five counts, reducing sensitivity to isolated high-energy outliers.

Reviewer comment:

“l. 288. The reference to electron-positron pair production is puzzling. It is unclear why this interaction is singled out, whereas photoelectric absorption and Compton scattering [...] are not mentioned.”

Response:

The passage has been revised to avoid undue emphasis on pair production and to provide a more balanced description of the relevant photon-interaction processes during post-field transport, including Compton scattering and photoelectric absorption.

The careful assessment of the manuscript and the constructive comments provided by Dr. Klimov are greatly appreciated.

Comment:

“For example, how significant is the discrepancy between the K coefficients given for different stations on page 6? It's clear they differ, but the calculation accuracy and method errors aren't provided. And if they're on the order of 10%, then the first two numbers (10.2 and 11.2) are barely distinguishable. So, it is better to demonstrate numbers with errors, as it is done in the second part of the work, where free path distance is discussed.”

Response:

The uncertainties of the station-specific K coefficients have now been explicitly evaluated and included in the revised manuscript. The effective avalanche length λ_eff is obtained from the slope of the linear fit to ln N_e(d), and the standard error of this slope is propagated first to λ_effand subsequently to the fitted coefficient K. The resulting values are

K_Aragats = (10.2 ± 0.4) × 10³ kV,

K_LHAASO = (11.2 ± 0.3) × 10³ kV,

K_{Lomnický Štít} = (8.9 ± 0.4) × 10³ kV,

K_{Nor Amberd} = (8.6 ± 0.2) × 10³ kV.

Some neighboring values show partial overlap within their fit-derived uncertainties. We therefore interpret the result as station-to-station variability within the atmospheric conditions considered, rather than claiming that every pair of station-specific coefficients is individually distinct. These uncertainties represent the statistical uncertainties from the fitting procedure and do not include possible systematic uncertainties associated with the atmospheric profile, field geometry, fitting interval, or CORSIKA physics settings.

Comment:

“A similar question concerns the value of the coefficient ‘a’, introduced on page 7. It turned out to be -0.92, but could it simply be -1? Which, in principle, might be logical: the normalized K coefficient should be divided by the density.”

Response:

The generalized empirical relation was fitted jointly to the combined dataset from all four stations, with K and a optimized simultaneously. The joint fit gives

K = (6.39 ± 0.26) × 10³ kV, a = −0.923 ± 0.083.

The fitted value of a is statistically consistent with a = −1. However, in the generalized relation

λ_eff= K · nᵃ / (E_z − E_th),

the density dependence of the runaway threshold is already included through E_th = E_₀ · n. The factor nᵃ therefore represents an additional empirical density dependence beyond the threshold scaling. Thus, although a = −0.923 ± 0.083 is consistent with −1, it should not be interpreted as a recovery of the classical inverse-density scaling. The fitted values of K and a are retained because they are obtained directly from the joint optimization of the simulation dataset.

Comment:

“Another point is: Is it possible to say something about correspondence of the spectra presented in fig 3 and 4 to the experimental ones, measured at the stations? It is beyond the scope of the paper, but it is interesting to know, how realistic these simulated spectra.”

Response:

A direct quantitative comparison with measured spectra would require information on the actual thunderstorm electric-field geometry, field termination height, detector response, and event-specific atmospheric conditions. These quantities are not constrained within the present simulation setup, and such a comparison is therefore beyond the scope of the present study.

The spectra referred to as Figs. 3 and 4 in the original manuscript are presented as Figs. 5 and 6 in the revised manuscript. They are used to characterize the simulated transport behavior as the field-to-detector distance varies, rather than to reproduce a particular observed TGE spectrum. The simulations show rapid attenuation of the electron component and more gradual attenuation of the gamma-ray component with increasing field-free propagation distance, providing the basis for the subsequent FPD analysis.

A quantitative simulation-to-observation comparison, including detector response and event-specific electric-field configurations, would require a dedicated study.

Comment:

“One more minor comment: lines 53–55 in the sentence ‘While this equation has been widely applied, it does not explicitly account for variations in AEF and atmospheric density and is often limited to narrow altitude and AEF ranges.’ it is better to give values of altitude ranges and AEF to which the equation is applied and give references.”

Response:

The corresponding discussion has been revised to avoid assigning a single altitude or electric-field range to the classical relation, because previous parameterizations were obtained under different electric-field and atmospheric-density conditions.

The revised Introduction now states that Dwyer (2003) reported K ≈ 7.2 × 10³ kV, while Coleman and Dwyer (2006) reported a related parameterization for electric fields above 300 kV/m and a different empirical fit for fields closer to the runaway threshold. Their scaling of the electric field and avalanche length to equivalent sea-level values is also discussed. The applicability of the avalanche-length parameterization is therefore described in relation to the electric-field range and atmospheric conditions for which the corresponding parameterization was obtained, with the relevant references included.