© R. Kataoka et al., Published by EDP Sciences 2022

1 Introduction

Measuring Galactic Cosmic Rays (GCR) is important, playing a unique role in diagnose the space weather and space climate. For example, the real-time data of worldwide neutron monitor network is utilized for predicting the radiation dose of aircrews during Ground-Level Enhancement (GLE) events (Kataoka et al., 2014; Sato et al., 2018). The long-term GCR data is also used to examine the solar cycle prediction of the aircrew dose (Miyake et al., 2017). Furthermore, for transient space weather events such as Coronal Mass Ejections (CMEs) passing through the Earth, we can also estimate the large-scale magnetic field structure of CMEs by analyzing the Global Muon Detector Network (GMDN) data (Kihara et al., 2021).

We started neutron and muon measurements at Syowa Station in the Antarctic (69.01° S, 39.59° E; the vertical cutoff rigidity is 0.4 GV) in February 2018 (Kato et al., 2021). The simultaneous measurements of neutron and muon flux offer a unique perspective relative to the other neutron monitors and muon detectors that do not share the same location. The median primary rigidities for the neutron monitor and muon detector (vertical) are 16.3 GV and 53.6 GV, respectively (Appendix A). We estimated the median primary rigidities by integrating the response functions of secondary neutrons and muons to primary cosmic rays (Nagashima et al., 1989; Murakami et al., 1979).

During the 4-year cosmic ray observation at Syowa Station, the solar activity gradually changes across the solar minimum. Figure 1 shows the solar minimum at the end of 24th and beginning of 25th 11-year solar cycles. The top panel shows the Sun’s polar magnetic field (http://wso.stanford.edu/). The north-south average of the magnetic field strength (green dots) started decreasing in late 2019. The second panel shows the sunspot number (https://www.sidc.be/silso/), passing through the solar minimum in late 2019. The bottom panel shows the neutron monitor count rate at Oulu (https://cosmicrays.oulu.fi/), starting to decrease in early 2020.

Figure 1 (a) The 10-day values of 30-day running averaged Sun’s polar magnetic field, (b) daily mean sunspot number, and (c) daily mean neutron monitor count rate at Oulu, Finland. The polar magnetic field strength (green dots) started decreasing in late 2019, while the neutron count rate started decreasing in early 2020.

This paper aims to introduce new methods for corrections of temperatures and snow cover effects using the first 4-year data of neutron count rate and muon count rate at Syowa Station. Specifically, the temperature correction is essential for muons, while the snow cover correction is vital for neutrons (Bütikofer, 2017). Section 2.2 proposes a new method for the snow cover correction based on the analytical radiation model. In Section 3, we show comparison of machine learning approach to other methods that are more rooted in physics of the processes involved and we discuss that a machine-learning technique combined with meteorological reanalysis data can be useful for these different corrections for muons and neutrons. Finally, concluding remarks are summarized in Section 4.

2 Methods of analysis

2.1 Mass weighted temperature correction for muon count rate

The muons are created in the air shower developing in the atmosphere. The atmospheric density and temperature profiles affect the production and propagation of muons. Due to the short lifetime of muons, muon count rate at ground level depends on the distance between production layer and detection plane. The count rate decreases with the increase of this distance due to atmospheric expansion caused by the increase of atmospheric temperature. Therefore, applying the so-called temperature corrections to muon data is necessary before analyzing the primary cosmic ray variations.

The Mass Weighted (MSS) method (Mendonsa et al., 2016) first calculates the mass-weighted temperature T_MSS from the temperatures T_i at an atmospheric layer at altitude h_i as

$$T_{\mathrm{MSS}}=\sum _{i=0}^{n}w\left(h_i\right)T\left(h_i\right).$$(1)

The air-mass weight function w is defined as

$$w\left(h_i\right)=(x\left(h_i\right)-x(h_{i+1}\left)\right)/x\left(h_0\right),$$(2)where x is the atmospheric depth in unit of (g/cm²) and h₀ is the ground-level altitude. In this paper we use the atmospheric pressure at h_i for x(h_i) in equation (2).

The correction based on the linear correlation between the T_MSS and muon count rate is the best method for muon detector data (Mendonsa et al., 2016). In this study, we applied the same MSS method to the muon detector data at Syowa Station as Kato et al. (2021) documented. To calculate the T_MSS, the temperature data was obtained from GDAS (Global Data Assimilation System; https://www.ncei.noaa.gov/products/weather-climate-models/global-data-assimilation) at 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, and 20 hPa. Kato et al. (2021) reported how the coefficient for MMS is calculated with the obtained values. The corrections of barometric effects for both muons and neutrons are also reported by Kato et al. (2021).

2.2 PARMA-based snow cover correction for neutron counting rate

The pressure correction is the major atmospheric correction necessary for the neutron count rate because the lifetime of neutrons is long and the air-shower production dominantly depends on the air mass above the detector. It is, however, also known that the neutron monitor count rate decreases when the accumulation of snow increases around the detector housing because of albedo neutrons scattered from the soil moisture or the snow cover of the surrounding ground area (Schron et al., 2016; Brall et al., 2021). The snow cover effect is not negligible for the NM64-type detector used at Syowa Station because the thickness of the reflector is only 7.5 cm (c.f., 28 cm in the IGY detector), and the evaporation neutrons produced from the surrounding material can significantly contribute to the counting rate with ~5% (Hatton, 1971). This paper discusses the snow cover effect in the polar region, while Ruffolo et al. (2016) discussed a similar “water vapor” environmental effect in the tropical region. Here, we do not anticipate the impact of the snow accumulation on the roof because the strong wind at Syowa Station tends to blow the roof snow away in a short time scale.

Japan Meteorological Agency observes the snow cover depth at Syowa Station at the Kitano-Ura area, several hundred meters away from the cosmic ray detector. The actual snow cover around the neutron monitor can be different depending on neighbor buildings and ground slopes. Also, the most significant human activity in snow removal is usually in November or December, depending on the occurrence of blizzard activities, while minor removals have been done randomly throughout the year. Therefore, to estimate the snow cover depth around the detector, we first use the observed snow cover depth between February and November and reset the estimated snow cover depth to be zero at the beginning of February each year. We then linearly interpolated the snow cover depth between November and February. Further, we put zero in the data gap of the snow cover estimation before June 2018, considering the microgravity observation at that time (Aoyama et al., 2016).

For the snow cover correction, we adopted the PARMA (PHITS-based Analytical Radiation Model in the Atmosphere) model (Sato, 2015, 2016) for calculating the neutron count rates. It can reproduce the influence of the surrounding environment on the neutron flux as a function of the underground water density (Sato & Niita, 2006). The calculated neutron fluxes were converted to the count rates, using the response function of the NM64-type detector (Sato et al., 2014). In this study, we assumed the linear relationship between the estimated snow depth, d, and the underground water density supplied to the PARMA model, q, as written by q = c₁ + c₂d, where c₁ and c₂ are the constant parameters. We fixed the numerical value of c₁ to be 0.20, commonly used for calculating the neutron monitor count rate (Sato, 2015). On the other hand, we regarded c₂ as a free parameter determined to minimize the χ² value between the calculated and measured count rates. We can estimate the snow correction factor from the ratio between neutron count rates calculated by PARMA with and without the snow cover effect (Appendix C).

Figure 2 shows the observed snow cover depth, the pressure-corrected relative count rates of the neutron monitor obtained from the observation, and the results of PARMA with different c₂ parameters. The relative count rates are normalized to their mean value in February 2019, when there was little snow around the neutron monitor. It is apparent from Figure 2 that the measured count rates are anti-correlated with the snow cover depth. The best-fit value of c₂ to the observation is 0.45 m⁻¹, with the maximum coefficient of determination R² = 0.69.

Figure 2 Observed snow cover depth (top panel) and the pressure-corrected relative count rates of the neutron monitor as obtained from the observation and the PARMA calculation with different c2 parameters (bottom panel).

2.3 ESN for muon and neutron counting rates

The Echo State Network (ESN) is a kind of recurrent neural networks (Jaeger, 2001), known as a novel method that works with a relatively small set of training data, which is especially suitable to the analysis of times series. The basic structure of the ESN model used in this paper is essentially the same as was used by Kataoka and Nakano (2021). Therefore, only the essential part is repeated below. The input vector u, reservoir state vector x, and output vector y are defined by N-points time series of n = 1, 2, 3…N as follows:

$$\mathit{u}\left(n\right)=\left(\begin{array}{c}{u}_1\left(n\right)\\ ⋮\\ u_{\mathrm{Nu}}\left(n\right)\end{array}\right),\hspace{0.5em}\mathit{x}\left(n\right)=\left(\begin{array}{c}{x}_1\left(n\right)\\ ⋮\\ x_{\mathrm{Nx}}\left(n\right)\end{array}\right),\hspace{0.5em}\mathit{y}\left(n\right)=\left(\begin{array}{c}{y}_1\left(n\right)\\ ⋮\\ y_{\mathrm{Ny}}\left(n\right)\end{array}\right).$$(3)

The reservoir state vector x consists of a large number of nodes, which are updated in time with the input vector u and the previous state of the nodes as follows:

$$\mathbf{x}\left(n+1\right)=f\left(W^{\mathrm{in}}\mathbf{u}\left(n+1\right)+W\mathbf{x}\left(n\right)\right) \left(n=0, 1, 2,\dots \right),$$(4)where we use hyperbolic tangent as the function f and fix the weight matrices Wⁱⁿ and W. To make the random and sparse node connections of W, we set the number of nodes, Nx, to be 10³, where only 10% of the matrix elements are random values between −1.0 and 1.0, and the rest 90% are zero. We selected the spectral radius (maximum eigenvalue) of W below unity to satisfy the echo state property that guarantees the independence of the reservoir state to the initial values (Jaeger, 2001). In this paper, we optimized the spectral radius to be 0.95 (Appendix B). The output vector y is calculated by the linear combination of the output weight matrix and the reservoir state vector as follows:

$$\mathbf{y}\left(n+1\right)=W^{\mathrm{out}}\mathbf{x}\left(n+1\right)\hspace{0.5em}\left(n=0, 1, 2,\dots \right),$$(5)where W^out is the output weight matrix. We train only the output weight matrix W^out by the set of T-point time series of input vectors X and desired output vectors D:

$$X=\left[\mathbf{x}\left(1\right),\dots ,\mathbf{x}\left(T\right)\right],\hspace{0.5em}D=\left[\mathbf{d}\left(1\right),\dots ,\mathbf{d}\left(T\right)\right].$$(6)

The least-squares method to minimize the difference between the outputs y and d can be represented by a standard linear regression as follows:

$$W^{\mathrm{out}}=D{X}^T{\left(X{X}^T\right)}^{-1}.$$(7)

The temperature data was obtained from GDAS at 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, and 20 hPa. In the above ESN model (Kataoka & Nakano, 2021), we use temperature data as the 15-dimensional input vector (Nu = 15) and count rate as the one-dimensional output vector (Ny = 1). In this study, about four-year data (from February 1, 2018 to March 31, 2022) are used for machine learning. We used six-hour averaged values, i.e., a total of 6080 data points exist in the time axis. The results using 24-hour values (a total of only 1520 data points) are also shown in Appendix (Figs. B2 and B3).

Usually, for many machine learning applications, we split the input dataset into training and testing data sets. In this paper, however, we used a different approach to demonstrate the possible method with limited dataset. We simply used the whole data for training only. The constructed model is then used for reproducing the training data, and not used for predicting any testing data. For longer term dataset with somewhat different purposes, however, it is possible to separate the training and testing datasets to examine the metrics of the constructed model. That can be a good future work.

3 Comparison among different correction methods

We compare the MSS and ESN temperature correction methods on the muon count rate (Fig. 3). Figure 3b shows the temperature effects on the muon count rate reproduced by the two methods. Note that a data gap in muon count rate existed at the beginning of 2021. Figure 3c shows the corrected count rates estimated from the ESN and MSS methods. The MMS corrected muon data is obtained from Kato et al. (2021). It is apparent that both methods similarly remove the seasonal variation due to the temperature effect. Similar decreasing trends of ~2% per year appeared in 2020 and continued in 2021 in Figure 3c. Although the seasonal variation remains in the corrected count rates, the start timing of the decreasing trend roughly corresponds to the solar minimum in late 2019.

Figure 3 Data-model comparison at Syowa Station. (a) Temperatures at 925, 700, 500, 400, 250, 150, 70, 50, and 20 hPa, (b) muon count rate in the unit of counts/min, (c) corrected count rate.

There are discrepancies between modeled and measured data as seen in Figure 3b: Namely, in 2019 both methods predict lower count rates than observation, while in 2021 and 2022 both methods overestimate the actual muon count rate. This is an expected result from the initial deterioration of muon detector over yearly time period. Such a discrepancy can be overcome if the corrections are calculated for shorter periods. However, in this study aiming at analyzing the yearly change of the mean count rate, we simply use the whole data as the first step.

For the snow cover correction of the neutron counts at Syowa Station, the same ESN method may also work for the neutron monitor at Syowa Station, because the time variation of the atmospheric parameters may control the resultant snow cover around the neutron monitor. However, the chain of physics and related human activities of snow removal is very complex.

Figure 4 shows the snow cover corrections applied to the neutron count rate at Syowa Station, using the PARMA and ESN models. The snow cover effects reproduced by two methods are shown in Figure 4b, while Figure 4c shows the corrected count rates (Appendix C). It is apparent that both methods similarly remove the seasonal variation due to the snow cover. In Figure 4c, a negative spike in May 2019 is a natural variation associated with multiple CMEs passing across the Earth. Similar Forbush decrease events associated with CMEs can be identified also in November 2021 and March 2022. On the other hand, the large bipolar variation in December 2020 is associated with the human activity of snow removal, which is not a natural variation of primary cosmic rays. Note also that the solar activity was relatively quiet in December 2020.

Figure 4 Data-model comparison at Syowa Station. (a) Temperatures at 925, 700, 500, 400, 250, 150, 70, 50, and 20 hPa, identical to Figure 3a, (b) neutron count rate in the unit of counts/min, (c) corrected count rate.

Figure A1 (a) Differential and (b) integral response functions of the muon detector (solid blue circles) and neutron monitor (solid orange diamonds) at Syowa Station. Differential integral response functions are normalized.

The decreasing trend starts in the somewhere early half of 2020 in both PARMA and ESN results, which is roughly consistent with Oulu neutron monitor data (Fig. 1c). We must monitor the actual snow cover situation around the neutron monitor in future observations for further detailed snow cover correction.

Note that the ESN-based correction of the snow-cover effect on the neutron count rate has a complex meaning. First, the input temperature data of GDAS correlates with the snow cover, which can also be confirmed by the ESN method (Appendix D). The agreement between the modeled and observed snow cover data is not very high, especially from the beginning of year 2020 onward. Therefore, in future work, some additional strong predictors would be needed as additional input variables, to increase the reliability of this method and more fully exploit the advantages of machine learning approach. Second, snow-drift accumulation effects, blizzard occurrence, and artificial snow removal activities are also related to the temperature data. The ESN learned the whole correlations among all of these natural and artificial effects.

If the snow data is available, we can use the PARMA-based correction for other similar neutron monitors. If not, we can still apply the ESN method using meteorological dataset. The results can then be compared among multi-stations in polar regions to examine the best possible corrections.

The decreasing trend started earlier in the muon data around the solar minimum in late 2019, while it started several months later in the neutron monitor data in early 2020 (see Figs. 3c and 4c). It has been established by observations that the temporal variation of the GCR intensity lags behind the solar parameters, such as the sunspot number, the interplanetary magnetic field magnitude at Earth, the neutral sheet tilt angle and the open solar magnetic flux (e.g., Cane et al., 1999; Koldobskiy et al., 2022), partly due to the GCR propagation time in the heliosphere which is shorter in high energy GCRs than in low energy GCRs. Recent numerical study demonstrated that the propagation time also depends on particles’ charge and the drift cycle (Strauss et al., 2012). However, the calculated propagation time is below 10 days for 10 GeV protons, much shorter than the lag seen in Figures 3 and 4.

By analyzing the long-term variations of neutron and muon count rates, Nagashima and Morishita (1980) reported that the time lag of GCR intensity variation behind the sunspot number is as long as 10 months and is systematically longer in odd solar activity cycles, including A > 0 to A < 0 transition of the solar polar magnetic field polarity than in even cycles including A < 0 to A > 0 transitions. Such long time lag and its 22-year variation are both verified by a recent work (Koldobskiy et al., 2022). The present paper documented an example only in A > 0 solar minimum. We have to examine the modulation in the A < 0 period and in different solar activity phases to investigate the energy-dependent solar modulation in more detail. As a future work, we can better address the hysteresis effects and the energy-dependent solar modulation by carefully applying similar correction methods as developed in this study to long-term data of both neutron monitors and muon detectors.

4 Conclusions

We proposed a new snow cover correction method for the neutron count rate using the PARMA model. We then showed that the ESN model combined with the GDAS temperature time series could be useful for the snow cover correction of neutrons and the temperature correction for muons, showing the reasonable agreement among different correction methods. From the comparisons of the corrected count rates, we conclude that muons likely started to decrease at least a few months earlier (late 2019) than neutrons (early 2020) following the onset of 25th 11-year solar cycle in 2019, which can be interpreted by a standard understanding of the energy-dependent intrusion of cosmic ray protons.

Acknowledgments

RK thanks Okinawa Institute of Science and Technology for hosting his sabbatical visit. RK thanks Gen Hashida, Shuki Ushio, Yuichi Aoyama, and Akihisa Hattori, for their advice on using the snow cover data at Syowa Station. We thank Japan Meteorological Agency for providing snow cover depth data at Syowa Station at the Kitano-Ura area. We acknowledge NOAA/NCEP for providing GDAS data (ftp://ftp.arl.noaa.gov/archives/reanalysis/). The polar magnetic field data is obtained from Wilcox Solar Observatory. Sunspot number is obtained from SILSO ver 2.0 (https://www.sidc.be/silso/datafiles). We thank the University of Oulu for providing neutron monitor data (https://cosmicrays.oulu.fi/readme.html). The editor thanks Du Toit Strauss and an anonymous reviewer for their assistance in evaluating this paper.

Appendix A

Response function and primary median energy

We calculate the response function for the Syowa neutron monitor using the method given by Nagashima et al. (1989). For the muon detector, for the vertical incident direction, we used the interpolation values from the table of the response function calculated by Murakami et al. (1979). We define the median primary rigidities as the rigidity at which the normalized integral response function equals 0.5.

Appendix B

Hyperparameter survey for echo-state network and the time resolution dependence

Root mean square error between the observation and model decreases as the spectral radius is greater for both muon modeling (left panel of Fig. B1) and neutron modeling (right panel of Fig. B1). To guarantee the echo state property, we selected the spectral radius of 0.95 (below unity). The examples using different time resolution of 24 h (less data points) are shown in Figures B2 and B3.

Figure B1 Root mean square errors between observed and modeled muon count rate (left) and neutron count rate (right), changing the spectral radius of the Echo State Network. We set the number of nodes to be 103. Selected temporal resolutions are 6 h (N = 6080) and 24 h (N = 1520).

Figure B2 Same as Figure 3 except for the time resolution, using 24 h values (N = 1520).

Figure B3 Same as Figure 4 except for the time resolution, using 24 h values (N = 1520).

Figure D1 (a) Input vectors of GDAS temperatures at 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, and 20 hPa and (b) output vector of the snow cover depth with the ESN result (red). Black curve is the observed snow cover in Figure 2a.

Appendix C

Corrected count rate and correction factor

The corrected count rate I_corr and correction factor A_corr can be defined as follows:

$$I_{\mathrm{corr}}=I_{\mathrm{obs}}{A}_{\mathrm{corr}},$$

$$A_{\mathrm{corr}}=\frac{I_{\mathrm{base}}}{I_{\mathrm{actual}}},$$where I_obs is the observed count rate, I_base and I_actual are the model-estimated count rates for the base and actual conditions, respectively. For the snow cover correction using PARMA, I_base and I_actual were calculated by setting q = 0.2 and 0.2 + 0.45d, respectively. As another example, for the standard pressure correction of neutron count rate, A_corr can be described as follows:

$$A_{\mathrm{corr}}=\exp \left(-R_{\mathrm{model}}\right),$$

$$R_{\mathrm{model}}=\ln \left(\frac{I_{\mathrm{actual}}}{I_{\mathrm{base}}}\right)=\frac{\beta }{100}∆P,$$where ΔP is the difference between the actual and base pressures, and β is the correlation coefficient in the unit of [%/hPa].

Appendix D

Example of the ESN model to reproduce the snow cover

To test the ability of the ESN model to reproduce the snow cover, we replace the output vector of the ESN model with the estimated snow cover data. The input data is the same GDAS temperature time series used in the main contents. The basic seasonal variation of the estimated snow cover is reasonably reproduced, although the disagreement is considerable in 2021.

References

All Figures

	Figure 1 (a) The 10-day values of 30-day running averaged Sun’s polar magnetic field, (b) daily mean sunspot number, and (c) daily mean neutron monitor count rate at Oulu, Finland. The polar magnetic field strength (green dots) started decreasing in late 2019, while the neutron count rate started decreasing in early 2020.
In the text

	Figure 2 Observed snow cover depth (top panel) and the pressure-corrected relative count rates of the neutron monitor as obtained from the observation and the PARMA calculation with different c2 parameters (bottom panel).
In the text

	Figure 3 Data-model comparison at Syowa Station. (a) Temperatures at 925, 700, 500, 400, 250, 150, 70, 50, and 20 hPa, (b) muon count rate in the unit of counts/min, (c) corrected count rate.
In the text

	Figure 4 Data-model comparison at Syowa Station. (a) Temperatures at 925, 700, 500, 400, 250, 150, 70, 50, and 20 hPa, identical to Figure 3a, (b) neutron count rate in the unit of counts/min, (c) corrected count rate.
In the text

	Figure A1 (a) Differential and (b) integral response functions of the muon detector (solid blue circles) and neutron monitor (solid orange diamonds) at Syowa Station. Differential integral response functions are normalized.
In the text

	Figure B1 Root mean square errors between observed and modeled muon count rate (left) and neutron count rate (right), changing the spectral radius of the Echo State Network. We set the number of nodes to be 103. Selected temporal resolutions are 6 h (N = 6080) and 24 h (N = 1520).
In the text

	Figure B2 Same as Figure 3 except for the time resolution, using 24 h values (N = 1520).
In the text

	Figure B3 Same as Figure 4 except for the time resolution, using 24 h values (N = 1520).
In the text

	Figure D1 (a) Input vectors of GDAS temperatures at 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, and 20 hPa and (b) output vector of the snow cover depth with the ESN result (red). Black curve is the observed snow cover in Figure 2a.
In the text