Local environmental effects on cosmic ray observations at Syowa Station in the Antarctic: PARMA-based snow cover correction for neutrons and machine learning approach for neutrons and muons

– Solar modulation of galactic cosmic rays around the solar minimum in 2019 – 2020 looks different in the secondary neutrons and muons observed at the ground. To compare the solar modulation of primary cosmic rays in detail, we must remove the possible seasonal variations caused by the atmosphere and surrounding environment. As such surrounding environment effects, we evaluate the snow cover effect on neutron count rate and the atmospheric temperature effect on muon count rate, both simultaneously observed at Syowa Station in the Antarctic (69.01 (cid:1) S, 39.59 (cid:1) E). A machine learning technique, Echo State Network (ESN), is applied to estimate both effects hidden in the observed time series of the count rate. We show that the ESN with the input of GDAS data (temperature time series at 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, and 20 hPa) at the local position can be useful for both the temperature correction for muons and snow cover correction for neutrons. The corrected muon count rate starts decreasing in late 2019, preceding the corrected neutron count rate which starts decreasing in early 2020, possibly indicating the rigidity-dependent solar modulation in the heliosphere.


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.
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 The air-mass weight function w is defined as where x is the atmospheric depth in unit of (g/cm 2 ) and h 0 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-climatemodels/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).

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(Sato, , 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 1 + c 2 d, where c 1 and c 2 are the constant parameters.We fixed the numerical value of c 1 to be 0.20, commonly used for calculating the neutron monitor count rate (Sato, 2015).On the other hand, we regarded c 2 as a free parameter determined to minimize the v 2 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 pressurecorrected relative count rates of the neutron monitor obtained from the observation, and the results of PARMA with different c 2 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 2 to the observation is 0.45 m À1 , with the maximum coefficient of determination R 2 = 0.69.

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: 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: where we use hyperbolic tangent as the function f and fix the weight matrices W in and W. To make the random and sparse node connections of W, we set the number of nodes, Nx, to be 10 3 , 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: 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: The least-squares method to minimize the difference between the outputs y and d can be represented by a standard linear regression as follows: 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.

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 R. Kataoka et al.: J. Space Weather Space Clim. 2022, 12, 37 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.
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.
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 R. Kataoka et al.: J. Space Weather Space Clim. 2022, 12, 37 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, snowdrift 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 energydependent solar modulation in more detail.As a future work, we can better address the hysteresis effects and the energydependent solar modulation by carefully applying similar correction methods as developed in this study to long-term data of both neutron monitors and muon detectors.

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.
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: where DP is the difference between the actual and base pressures, and b is the correlation coefficient in the unit of [%/hPa].R. Kataoka et al.: J. Space Weather Space Clim. 2022, 12, 37

Figure 1 .
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.

Figure 2 .
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 c 2 parameters (bottom panel).

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