A method is developed to rank Forbush decreases (FDs) in the galactic cosmic ray radiation according to their expected impact on the ionization of the lower atmosphere. Then a Monte Carlo bootstrap‐based statistical test is formulated to estimate the significance of the apparent response in physical and microphysical cloud parameters to FDs. The test is subsequently applied to one ground‐based and three satellite‐based data sets. Responses (>95%) to FDs are found in the following parameters of the analyzed data sets. AERONET: Ångström exponent (cloud condensation nuclei changes), SSM/I: liquid water content, International Satellite Cloud Climate Project (ISCCP): total, high, and middle, IR‐detected clouds over the oceans, Moderate Resolution Imaging Spectroradiometer (MODIS): cloud effective emissivity, cloud optical thickness, liquid water, cloud fraction, liquid water path, and liquid cloud effective radius. Moreover, the responses in MODIS are found to correlate positively with the strength of the FDs, and the signs and magnitudes of the responses agree with model‐based expectations. The effect is mainly seen in liquid clouds. An impact through changes in UV‐driven photo chemistry is shown to be negligible and an impact via UV absorption in the stratosphere is found to have no effect on clouds. The total solar irradiance has a relative decrease in connection with FDs of the order of 10 −3 , which is too small to have a thermodynamic impact on timescales of a few days. The results demonstrate that there is a real influence of FDs on clouds probably through ions.

1 Introduction Galactic cosmic rays (GCRs) are thought to affect cloud cover on Earth, through ionization and subsequent effects on aerosol processes [Dickinson, 1975; Svensmark and Friis‐Christensen, 1997; Marsh and Svensmark, 2000; Bagó and Butler, 2000]. A prevalent approach to evaluating this idea has been to investigate effects of coronal mass ejections—events where a cloud of magnetized plasma is ejected from the solar corona and travels out into interplanetary space. This ejected plasma cloud tends to screen out galactic cosmic rays from its interior and if the plasma cloud hits Earth, it may result in a sudden decrease in the amount of cosmic rays reaching the atmosphere as measured by neutron monitors. Such an event, known as a Forbush decrease (FD) [Forbush, 1937], may evolve within hours, with the depression in cosmic ray counts recovering over a week or more, as the plasma cloud passes Earth and continues farther out from the Sun. An example of a FD is shown in Figure 1. Fewer cosmic rays generate less ionization in the atmosphere, such that these events present an opportunity for testing the cosmic ray/cloud link. Figure 1 Open in figure viewer PowerPoint The large Forbush decrease of 2003 with the minimum in the daily averaged signal occurring on 31 October as recorded by the CLIMAX neutron monitor. The black and red lines, respectively, correspond to the hourly and daily averaged values. The latter is used in this work. N R is the reference level before the FD, while N min =N(t 0 ) is the daily averaged minimum value during the FD (see text). Svensmark et al. [2009] found a significant signal in both aerosols (specifically the aerosol Ångström exponent) and liquid, low clouds. However, a debate is still ongoing with regard to whether there actually is an effect. Studies by Sloan and Wolfendale [2008], Laken et al.[2009], Calogovic et al. [2010], and Laken and Calogovic [2011] found no statistically significant signal during or following FDs. Kristjánsson et al. [2008] and Todd and Kniveton [2004] found some response in cloud satellite data, while Pudovkin and Veretenenko [1995], Harrison and Stephenson [2006], and Harrison and Ambaum [2010] found a response based on surface observations. A signal in middle‐ and high‐level clouds has been reported by Rohs et al. [2010], and Dragić et al. [2011] found a signal by using the diurnal temperature range as a cloud proxy, which allowed them to extend the data range farther back than the onset of satellite measurements. Experimentally, it has been shown that ions do promote the formation of new small (3 nm sized) aerosols [Svensmark et al., 2007; Kirkby et al., 2011], and one experiment suggests that ions also help the growth of aerosols to cloud condensation nuclei sizes (>50 nm) [Svensmark et al., 2013]. However, from a modeling point of view, it is uncertain whether a variation in ion‐induced nucleation may translate into an observable change in cloud condensation nuclei (CCN) and thus in clouds. Bondo et al. [2010] suggest that an aerosol effect could be observable under atmospheric conditions while general circulation modeling gives rise to much smaller responses in the CCNs [Kazil et al., 2006; Pierce and Adams, 2009; Snow‐Kropla et al., 2011; Yu and Luo, 2014]. Of these studies, Yu and Luo [2014] find a response almost an order of magnitude larger than previous estimates but still insufficient to explain the large observed variations in the ocean heat content over the solar cycle [Shaviv, 2008; Howard et al., 2015]. The above uncertainties make it is desirable to use observations to better constrain a possible effect. The present work aims to improve on previous statistical strategies, as well as apply these to both previously examined and unexamined cloud and aerosol data. We consider in depth the behavior of several cloud parameters during FDs and compare with cloud theory. We also improve on the treatment of potential effects from total solar irradiance (TSI) and ultra violet (UV) light during FDs by taking a solar spectral approach, rather than a UV tracer approach as done in previous studies. Four independent atmospheric data sets are used and described in section 2. They include a land‐based set of measurements of aerosol optical thickness probing the amount of small aerosols and three distinct satellite‐based cloud data sets measuring a range of physical and microphysical cloud parameters. Data from one or more of the above four independent sources have all been used in previous studies [Todd and Kniveton, 2004; Kristjánsson et al., 2008; Sloan and Wolfendale, 2008; Svensmark et al., 2009; Laken et al., 2009; Calogovic et al., 2010]. The response in atmospheric parameters is expected to depend on the strength of the FDs, namely, if a given FD produces a larger change in atmospheric ionization, it should correspondingly induce a relatively large response in the atmospheric variables. A method to estimate the strength of individual FDs was used implicitly in Svensmark et al. [2009]. In order to make the present paper self‐contained, the full method is presented in section 3. The statistical model is defined in section 4. This model is based on Monte Carlo bootstrap statistics, and it takes the finite growth time of aerosols into account. This statistical model is applied to each of the four data sets in section 5, and the resulting statistical significance is then determined. The rich Moderate Resolution Imaging Spectroradiometer (MODIS) data allow for the study of six relevant cloud parameters. If there exists a link between cosmic rays, aerosols, and clouds, a change in several cloud parameters can reasonably be expected on a global scale during a FD, where the sign and magnitude of the response should concur with expected cloud microphysics. In section 6, the expected parameter changes are calculated and compared to the results of the MODIS data set. Last, the results are discussed in section 7 and placed in context of possible mechanisms involving changes in the GCR flux, the Total Solar Irradiance (TSI), or UV light.

2 Data A number of data sets related to cosmic rays, total solar irradiance, ultraviolet radiation, atmospheric aerosols, and global cloud cover serving as the basis for the present analysis are presented below. 2.1 Neutron Monitor and Muon Data Cosmic ray variations at Earth have been monitored by either neutron monitors or by muon detectors. Both neutrons and muons are among the secondary particles produced when primary cosmic ray particles interact with atomic nuclei in the atmosphere, and the aforementioned detectors may therefore be used for probing the primary cosmic ray spectrum. In this study we employ neutron monitor data from all available stations within the temporal range 1987–2014. The number of stations amounts to about 130 with cutoff rigidity in the range 0–47 GV. Neutron monitor data can be obtained from the World Data Center for Cosmic Rays (WDCCR, http://center.stelab.nagoya-u.ac.jp/WDCCR/). Muon data are from the Multi‐Directional Cosmic Ray Muon Telescope at Nagoya and covers a cutoff rigidity of 60–119 GV. (http://www.stelab.nagoya-u.ac.jp/ste-www1/div3/muon/muon1.html). 2.2 Solar Spectrum Data For the analysis of the temporal variations in the solar electromagnetic spectrum we employ data from either composite solar spectral irradiance in the wavelength range 120–400 nm covering the time period 8 November 1978 to 1 August 2005 [DeLand and Cebula, 2008] or the Solar Radiation and Climate Experiment (SORCE) in the wavelength range 115–2416 nm covering the temporal period 14 April 2003 to 24 August 2015 (http://lasp.colorado.edu/home/sorce/data/). Furthermore, we use total solar irradiance (TSI) data from the VIRGO Experiment on the cooperative ESA/NASA Mission SOHO (version d41_61_0803) from PMOD/WRC. 2.3 Atmospheric Data observational data on aerosols in the atmosphere obtained from the solar photometers of the aerosol robotic network (AERONET) program [Schuster et al., 2006 observations of cloud liquid water content (CWC) over the world's oceans observed by the Special Sounder Microwave Imager (SSM/I) [Wentz, 1997 Weng et al., 1997 IR‐detected measurements of high, mid, low, and total IR clouds from the International Satellite Cloud Climate Project (ISCCP) [Rossow and Schiffer, 1991 daily observations of six key cloud parameters measured by the Moderate Resolution Imaging Spectroradiometer (MODIS) [Salomonson et al., 1989 ε), (2) the cloud optical thickness (τ), (3) the liquid water cloud fraction (CF, previously examined in Svensmark et al. [ 2009 R eff ), all from the “MOD08_D3” product. The ice cloud fraction is also used to a smaller degree. The atmospheric data used in this study are

3 Variations in the Primary Cosmic Ray Spectrum Caused by Forbush Decreases Measurements of the temporal evolution of the cosmic ray flux are primarily made by neutron monitors and to a lesser extent by muon telescopes. Since the response of neutron monitors depends on the location (latitude, longitude, and altitude) of the monitor, the monitors are sensitive to different parts of the primary cosmic ray energy spectrum. This fact will be used in the following to extract information on the variation in the primary cosmic ray energy spectrum during a FD. The primary spectral changes are important since they can be used to determine the changes in ionization throughout Earth's atmosphere, and thus be used to rank the strength of the FDs according to impact on ionization. 3.1 Extraction of the Primary Cosmic Ray Spectrum Individual FDs are identified over a period of almost three decades together with responses in the global network of about 130 neutron monitors (NMs) and muons from the Multi‐Directional Cosmic Ray Muon Telescope at Nagoya. The change in the primary cosmic ray spectrum at Earth (at 1 AU) during a FD cannot be known from a single neutron monitor, as the response is related to an integral over the cosmic ray spectrum and the response function of the neutron monitor. It is, however, possible to extract variations in the cosmic ray spectrum using multiple cosmic ray detectors as will be shown in the following. N(t), that a NM registers depend on the altitude and geomagnetic position of the instrument and is given by (1) P is the rigidity defined as P = pc/q, p and q, respectively, are the momentum and charge of the primary cosmic ray particle, and c is the speed of light. P c is the cutoff rigidity due to the geomagnetic field and S(h,P) is the yield or response function (the average number of counts in the NM, located at a height h above sea level, due to a primary cosmic ray particle of rigidity P) [Clem and Dorman, 2000 J(P,t) is the differential rigidity spectrum at 1 AU as a function of time t. Here J(P,t) is the unknown function whose variation during a FD we are aiming at estimating. A neutron monitor counts mainly the neutrons that are produced in the secondary shower events following the nuclear interactions of a primary cosmic ray particle with an atom high in Earth's atmosphere. Count rates,), that a NM registers depend on the altitude and geomagnetic position of the instrument and is given bywhereis the rigidity defined asand, respectively, are the momentum and charge of the primary cosmic ray particle, andis the speed of light.is the cutoff rigidity due to the geomagnetic field and) is the yield or response function (the average number of counts in the NM, located at a heightabove sea level, due to a primary cosmic ray particle of rigidity) [].) is the differential rigidity spectrum at 1 AU as a function of time. Here) is the unknown function whose variation during a FD we are aiming at estimating. P m as (2) From the above equation, one can define the median rigidityasi.e., the rigidity below which the NM registers 50% of its counts. The median rigidity characterizes a NM, since it depends on the location. One feature of the median rigidity is that it changes through the solar cycle. This is not a serious problem in the present study, since FDs are most frequent around solar maximum. We, therefore, elect to use the median rigidity at solar maximum for all NMs. The median rigidity of the NM data is based on vertical cutoff rigidity estimates and ranges from ≈10 GV (South Pole stations) to ≈47 GV (Ahmedabad, India). The multidirectional muon detector at Nagoya has 17 different viewing angles and represents 17 different paths of the muons through the atmosphere and, therefore, 17 different response functions. The median rigidity range of the Muon Telescope ranges from 60 GV to 119 GV. Together, the NMs and muon detector used in the present work cover the range from 10 GV to 119 GV in median rigidity. The analysis is based on daily averages. t 0 over which the minimum number of counts N(t 0 ) occurs is found. Then, data from all operational neutron monitors at that particular date are used to obtain reference levels, N R , prior to the minimum of the FD, defined as (3) (4) (5) Forbush decreases are here identified in the South Pole neutron monitor data as events having a relative decrease of at least 7%. This neutron monitor is used since it has the smallest cutoff rigidity of all monitors and thus has the largest response. When a FD is identified, the dayover which the minimum number of counts) occurs is found. Then, data from all operational neutron monitors at that particular date are used to obtain reference levels,, prior to the minimum of the FD, defined asi.e., a 15 day average of the neutron counts ending 1 day prior to the minimum. From this, the change in neutron counts is defined asand the relative change is defined as N R , ΔN, and N(t 0 ) are illustrated in Figure 1 for the strong FD event of 31 August 2003, as registered by the CLIMAX NM. It is possible to extract δN j,k for each of the operational NMs, where the index j identifies the NM, and index k the particular FD event. δN in the primary spectrum for a sample FD event. The data points correspond to measurements by neutron monitors operational at the particular date as function of median rigidity, evaluated at solar maximum. The group of points above 60 GV are from the Multi‐Directional Cosmic Ray Muon Telescope at Nagoya. One can then fit the data points with the functional form (6) A k and the exponent γ k are determined by a least squares fit, while the index k refers to the specific FD being investigated, and P is measured in gigavolts. The functional form of δn k given in equation Ahluwalia and Fikani, 2007 A = 229 and γ =− 0.87. Figure 2 depicts the extracted relative changein the primary spectrum for a sample FD event. The data points correspond to measurements by neutron monitors operational at the particular date as function of median rigidity, evaluated at solar maximum. The group of points above 60 GV are from the Multi‐Directional Cosmic Ray Muon Telescope at Nagoya. One can then fit the data points with the functional formwhere the amplitudeand the exponentare determined by a least squares fit, while the indexrefers to the specific FD being investigated, andis measured in gigavolts. The functional form ofgiven in equation 6 is regularly used as a good approximation to FD modulations of the energy spectra []. The blue line in Figure 2 is the result of a least squares fit of the above functional form to the strong FD event of 31 October 2003, resulting in the parameters= 229 and=− 0.87. Figure 2 Open in figure viewer PowerPoint P m of the stations, for the strong FD of 31 August 2003 (black star symbols). The blue curve shows the least squares fitted empirical function of equation A = 229 and γ =− 0.87. The blue line can also represent an exact theoretical relative change in the primary spectrum as indicated by equation P m , calculated for a standard NM‐64 neutron monitor response function and are a test of the approximation given by equation Relative changes of counts by neutron monitors and the Nagoya Muon Telescope plotted against the median rigidityof the stations, for the strong FD of 31 August 2003 (black star symbols). The blue curve shows the least squares fitted empirical function of equation 6 , with the parameters= 229 and=− 0.87. The blue line can also represent an exact theoretical relative change in the primary spectrumas indicated by equation 8 . Blue symbols are the relative changes in cosmic ray counts as a function of the median cutoff rigidity, calculated for a standard NM‐64 neutron monitor response function and are a test of the approximation given by equation 13 . The good agreement between the points and the curve demonstrates that the approximation is adequate and thus suggests that changes in the primary cosmic ray spectrum during FDs can be derived and quantified using the proposed method (see text). J 0 (P,t) and the perturbed spectrum J k (7) (8) Using equation 6 , the relation between the primary unperturbed differential rigidity spectrum) and the perturbed spectrumsuch that the change in the primary differential rigidity spectrum is (9) (10) (11) (12) It is now possible to relate between the responses in neutron monitors and the primary spectrum. The reference level derived from the detector counts is given bywhereand the minimum count during a FD is given by, J R (P) = J 0 (P,t 0 ) the relative change δN (see equation (13) By letting) =) the relative change(see equation 5 ) becomes To derive the second line of equation 13, the integration limits are changed using equation 2, while it is noted in the third line that S(h,P)J R (P) does not change sign and that the mean value theorem for integrals can be used to take out of the integration, where P 0 is in the interval P c ≤P 0 ≤P m . Finally P 0 is approximated by P m . Using the final form of equation 13 it is then possible to extract any spectral changes between two different sets of observations on any timescale relevant for GCR monitors. In particular, one can compare observations between solar maximum and solar minimum during solar cycles number 21–23, and characterize the spectral changes over the solar cycle with A SC =336 ± 46 and γ SC =−1.10 ± 0.04 [Svensmark et al., 2009]. This measure can then be used to relate the spectral changes during each FD to the spectral changes over a solar cycle. We continue by testing the above approximations numerically. With the aid of a response function S(P,h) for an NM‐64 neutron monitor [Clem and Dorman, 2000; Flückiger et al., 2008], one can calculate the relative change in the monitor counts δN for a given A k and γ k , as a function of the monitor's cutoff rigidity, P c , using equations 9 and 12. The relative change, δN calculated from the neutron monitor's response, can then be compared with the relative change in the differential rigidity spectrum, , (equation 8). The solid line in Figure 2 denotes , while the blue dots correspond to the numerically determined relative decrease δN, as a function of the numerically determined P m , (using equation 2). Figure 2 demonstrates that the approximation is satisfactory to within a few percent. The largest errors are for the lower median rigidities. Although the above test was not performed for a muon detector response function, there is no reason to suspect that our approximations are less valid in that case, as the response functions are of a similar form. With a threshold of 7% in the South Pole neutron monitor, a total of 26 events are found between 1989 and 2005. For each event, all available neutron monitors together with the Nagoya muon telescope are used to estimate the relative maximal decrease during the FD as a function of the monitors stated median rigidity P m . With the methods outlined above, it is then possible to extract the change in the primary cosmic ray spectrum for each FD. Table 1 summarizes the 26 largest FD events identified in this way. We proceed in the next section to rank the events according to their change in atmospheric ionization. Table 1. The 26 Strongest FD Events Over the 1987–2007 Period, Sorted by Strength Order Date Decrease (%) A ±δA γ 1 31/10/2003 119 229 10/9 −0.87 ± 0.02 2 13/6/1991 87 121 4/4 −0.74 ± 0.01 3 19/1/2005 83 273 16/15 −1.09 ± 0.02 4 13/9/2005 75 233 34/33 −1.07 ± 0.04 5 15/3/1989 70 93 14/12 −0.72 ± 0.06 6 16/7/2000 70 131 7/7 −0.86 ± 0.02 7 12/4/2001 64 153 12/11 −0.96 ± 0.03 8 29/10/1991 56 83 4/4 −0.76 ± 0.02 9 9/7/1991 54 84 4/4 −0.78 ± 0.02 10 29/11/1989 54 173 13/12 −1.08 ± 0.03 11 10/11/2004 53 95 8/8 −0.84 ± 0.04 12 26/9/2001 50 203 16/15 −1.18 ± 0.03 13 25/3/1991 48 82 15/13 −0.82 ± 0.07 14 17/7/2005 47 147 14/13 −1.07 ± 0.04 15 25/9/1998 45 123 45/33 −1.01 ± 0.14 16 27/7/2004 45 97 7/7 −0.91 ± 0.03 17 10/9/1992 44 206 46/38 −1.24 ± 0.09 18 31/5/2003 44 61 3/3 −0.74 ± 0.02 19 25/11/2001 39 75 15/13 −0.87 ± 0.08 20 15/5/2005 38 132 16/14 −1.12 ± 0.05 21 28/8/2001 37 152 15/14 −1.19 ± 0.04 22 27/8/1998 36 38 24/15 −0.63 ± 0.21 23 10/5/1992 35 50 6/5 −0.75 ± 0.05 24 27/2/1992 33 30 2/2 −0.57 ± 0.03 25 18/2/1999 33 38 3/3 −0.66 ± 0.03 26 2/5/1998 28 55 6/5 −0.88 ± 0.04 3.2 Variation in the Atmospheric Ionization Heck et al., 1998 T and an incident angle from the zenith in the interval 0≤α≤70° are the initial conditions for the cascade. Note that Usoskin and Kovaltsov [ 2006 I(P,h), the average ionization energy deposited at various heights in the atmosphere, is obtained. For energies in the range of 1 GeV–1000 GeV, the ion production in the atmosphere becomes (14) I(P,h) is the ion pair production at height h, caused by a primary particle with rigidity P. J(P) is the differential cosmic ray rigidity spectrum, which here follows the near–solar minimum observations of the Bess spectrometer [Sanuki et al., 2000 (15) The (maximal) change in the differential rigidity spectrum during a FD can be extracted by applying the methods outlined in the previous section. With it, one can calculate the resulting change in the ion production throughout the atmosphere. This is carried out through a Monte Carlo simulation of primary CRs incoming at different energies and the resulting shower structure of secondary particles. Specifically, the evolution of the shower is calculated using the CORSIKA code [], where a primary cosmic ray proton of kinetic energyand an incident angle from the zenith in the interval 0≤≤70° are the initial conditions for the cascade. Note that] have developed a model of atmospheric ionization based on similar methods where they use zenith angles up to 90°. For any primary particle energy, 10,000 showers are calculated, and), the average ionization energy deposited at various heights in the atmosphere, is obtained. For energies in the range of 1 GeV–1000 GeV, the ion production in the atmosphere becomeswhere) is the ion pair production at height, caused by a primary particle with rigidity) is the differential cosmic ray rigidity spectrum, which here follows the near–solar minimum observations of the Bess spectrometer []. Thus, any change in ionization during FDs is probed relative to the spectrum close to solar minimum. The change in ionization due to a FD is then given by Figure 3 (top) shows the ion production as a function of altitude (see, e.g., Bazilevskaya et al. [2008] for comparison with observations). The black thick curve is the ion production at solar minimum estimated by the reference spectrum at solar minimum [Sanuki et al., 2000]. The thin curves are the reduction in ion production relative to the solar minimum spectrum due to each of the 26 FDs based on the fitted A and γ, as given in Table 1. The lowest thin red curve is the exceptional event of October 2003. Using the functional form of equation 6 describing the solar cycle modulation of cosmic rays, one obtains the thick red curve in Figure 3 for the spectrum at solar maximum. Figure 3 (bottom) contains the same information except that the reduction in ion production is normalized here to the reduction from solar maximum to solar minimum, i.e., the difference between the black curve and the red curve at each height h. The ion production is here calculated assuming a cutoff rigidity of 5 GV. Using different values will change the detailed shape of the individual FD curves. Due to the uncertainties in the A and γ parameters, there is also an uncertainty in the derived atmospheric ionization, typically ±5% in the y axis values in Figure 3 (bottom). This can cause the rankings of the FDs which are close in strength to shift. An improvement that could be made to the model would be to run it for other parts of the globe than the 45° latitude and 5 GV cutoff rigidity we use with the U.S. Standard Atmosphere. This would require using both atmospheres and cutoff rigidities for several latitudes. In section 5.4 we look at the connection between FD strength and response, which gives an indication on whether the rankings are reasonable. Figure 3 Open in figure viewer PowerPoint Ion production in the atmosphere. (top) The absolute ion production in a U.S. standard atmosphere as a function of altitude. The thick black curve is the ion production under solar minimum conditions and the thick red curve is during solar maximum, corresponding to latitude of approximately 45° and a cutoff rigidity of 5 GV. The individual thin lines represent the depression relative to the conditions of solar minimum due to the FD events given in Table 1 . The lowermost thin red line corresponds to the very strong FD event in October 2003. (bottom) The individual FDs normalized to the solar cycle variation. Solar minimum is given by the thick dashed black line and solar maximum by the thick dashed red line, which in Figure 3 (top) corresponds to the difference between the thick black red curves. Of special interest in this paper are the changes in the lowermost 3 km of the atmosphere here shown as the grey area.

4 Aerosol and Cloud Variations During Forbush Decreases The main goal in the present work is to study whether FDs have an impact on atmospheric aerosols and clouds. As already mentioned in section 1, a number of studies have looked into this question with diverging conclusions. It is therefore imperative to take a robust statistical approach. In this section such an approach is described. 4.1 Processing Forbush Decrease Time Series U(t), where t enumerates the day. Using all the FDs in Table U(t) as follows: A FD minimum date t 0 is obtained for FD number i of Table 1. A temporal range encompassing t 0 is set to 36 days, 15 days prior to t 0 , and 20 days after, i.e., U(t) with t = [t 0 −15,t 0 +20] days. Any linear trend and sample mean is removed over the 36 days. The result is referred to as one FD unit. The data under consideration in this study are discrete daily measurements of aerosol or cloud parameters covering a number of years. Any such data series will be denoted), whereenumerates the day. Using all the FDs in Table 1 , a collection of FD time series are compiled from) as follows: The collection of FD units is then denoted F i (t). The next step is to reduce the daily data in the days following the FD into a quantifiable measure of a possible response in the different time series. Two such measures are defined below and then also used as a basis for statistical tests. The growth time of small aerosols into cloud condensation nuclei (CCN) is expected from aerosol dynamics to be within a range of a few days to a maximum of a little more than a week, depending on the rate of aerosol growth [Kulmala et al., 2004]. Thus, when averaging clouds or aerosols over a large area of the Earth, a response in any measured parameter should be correlated in the days following the onset of the FD. This further suggests that an integral over a period of days following the FD should be a good measure of the response. Such a measurement would automatically deal with the autocorrelations in the time series. (16) w i are statistical weights with the property that , while N FD is the total number of FDs (depending on the temporal range of the particular data set). Three weight distributions are used (17) The first type of measure is a weighted sum of the FD units, defined aswhereare statistical weights with the property that, whileis the total number of FDs (depending on the temporal range of the particular data set). Three weight distributions are used (18) (19) The first distribution provides weight according to the derived strength of the FD (Table 1) and includes all FDs. The second distribution provides equal weight to the five strongest FDs and ignores all other. Finally, the third weight distribution gives equal weight to all FDs, strong or weak. Additional weight distributions can of course be applied, but they would not offer notably different information on the significance of any signal which may be found, since the above distribution covers the most extreme range of distribution, from just the strongest to all at equal weights. t 1 to day t 2 following the FD minimum, which is the temporal range over which a signal is expected to appear (20) The second measurement type is the response defined as the summed signal within the period from dayto dayfollowing the FD minimum, which is the temporal range over which a signal is expected to appear FS is therefore a number. The next step is then to define statistics for the significance of these measurements. 4.2 Bootstrap Samples and Statistical Measures It is always preferable to minimize the amount of assumptions concerning the nature of any system under examination, as it may end up biasing the result. In situations where the underlying statistical distribution of a variable is simply unknown, the bootstrap method serves as an excellent way to perform distribution‐independent statistics. Furthermore, the bootstrap method can also handle issues inherent to time series analysis, such as autocorrelation [Efron and Tibshirani, 1994]. U(t) as follows: A random date t 0 is drawn with replacement. A temporal range encompassing t 0 is set to 36 days, 15 days prior to t 0 , and 20 days after, i.e., U(t), where t = [t 0 −15,t 0 +20] days. Any linear trend and sample mean is removed over the 36 days. The result is referred to as one sample unit. Steps 1–3 are repeated until a number of sample units corresponding to the number of Forbush decreases (N FD ) has been collected. This collection is referred to as a bootstrap sample. The procedure described in steps 1–4 is repeated a number of times, based on the size of the data set being used. The final number of bootstrap samples is N B , and denoted B i,j (t), where i denotes the sample unit number, and j is the bootstrap sample number. The statistical samples are produced and processed in much the same way as the FD time series of section 4.1 ; however, using random samples from) as follows: N FD sample units, it is easy to apply the weights of equations B i,j (t), constructs similar to those of section (21) (22) j still indicates the bootstrap sample number. For more information on these constructs, see the previous section. Note that since each bootstrap sample containssample units, it is easy to apply the weights of equations (17)–(19)-(17)–(19) to the bootstrap samples as well. Using the randomly generated), constructs similar to those of section 4.1 can be produced:andwhere indexstill indicates the bootstrap sample number. For more information on these constructs, see the previous section. measured value time series and the expectation values of theN B bootstrap samples (23) (24) In the general case it is possible to define the test statistic as the difference between theand theand for the integrated signal, t) and XFS are drawn from the same distribution as the bootstrap samples, defined as (25) (26) In our specific case the two sums (the second terms) in equations 23 and 24 are zero by construction. The question now is whether the above observed values XFW() and XFS are drawn from the same distribution as the bootstrap samples, defined asand for the integrated signal, Efron and Tibshirani, 1994 XFW(t) is obtained for each time step over the 36 day period as, (27) This can be answered probabilistically by calculating an achieved significance level (ASL) [, chapter 15] from the bootstrap samples. The achieved significance for the time series) is obtained for each time step over the 36 day period as, t indicates a particular day in the time series and n the index of the bootstrap sample. In the case of the integrated response, XFS, the achieved significance level becomes (28) Whereindicates a particular day in the time series andthe index of the bootstrap sample. In the case of the integrated response,, the achieved significance level becomes Such is the basis of the model and FW(t) and FS can be readily calculated and their values compared with the empirically determined properties of the distribution function, from identical measurements, BW j (t) and BS j , performed on the bootstrap samples.

5 Statistical Test of Atmospheric Data During Forbush Decreases Here the statistical procedure outlined in the previous section will be applied to the four data sets. 5.1 Forbush Decreases and AERONET Data The AERONET observations from more than 700 stations are available over the period 1998–present, which spans 17 FDs from the list in Table 1. AERONET data provide information on the transmission of solar radiation from the top of the atmosphere down to the surface at a number of different wavelengths. To quantify the relative blocking of sunlight at different wavelengths, we use here the Angstrom exponent α in the aerosol extinction law, defined through τ(λ i ) = τ 1 λ−α, where τ(λ i ) is the aerosol optical thickness at the wavelength λ i and τ 1 is the approximate optical thickness at λ= 1 μm. The AERONET data give the fitted exponent α 1,2 at two wavelengths λ 1 and λ 2 which provide information about the relative abundance of fine aerosols. Long wavelengths respond to their volume fraction, while short wavelengths are more sensitive to the effective radius of the fine mode (<250 nm) aerosol [Schuster et al., 2006]. This motivates us in this study to use the shortest wavelengths at 340 and 440 nm based on the idea that if FDs have an effect on the aerosol production by decreasing nucleation, it should manifest itself as an increase in the effective radius of the aerosols. The measured time series of the Angstrom exponent for the wavelengths 340–440 nm (equation 23 using the weights from equation 18) is displayed in Figure 4 (top), for the five strongest FDs superposed with even weights. The effect of using different weights is explored in the following sections. It is seen that there is a reduction in the Angstrom exponent in the days following the FD events. The dotted lines denote the one and two sigma deviations (per day) calculated from the bootstrap samples. The hatched area is the period used to integrate the response (day 0 to day 8, where day 0 is the date of the FD minimum stated in Table 1). This period is chosen based on the period during which aerosols can grow to sizes detectable by AERONET, and is expected to be within a week [Kulmala et al., 2004]. Slight variations in the averaging periods used in the paper do not change the conclusions of the paper (see section 7.1). One should note that the AERONET measurements are all performed on land, and here aerosols grow faster than over the remote ocean due to higher condensable vapor pressures over land. Figure 4 Open in figure viewer PowerPoint AERONET response to FDs. (top) The superposed FD response in the AERONET data for the five strongest Forbush decreases (after 1998, when the data set starts), 15 days before the minimum in the cosmic ray flux and 20 days after. The dotted lines are one sigma and two sigma obtained using the MC analysis. The hatched area (day 0 to day 8) is the temporal interval used to make the integrated response. (bottom) The integrated interval of each of the bootstrap samples as a black dot. The red horizontal line is the size of the Forbush signal. The response in the integrated signal is at the 98.97% significance level. Figure 4 (bottom) shows the integrated values of the bootstrap time series samples as black dot symbols (equation 26). The red line shows the size of the signal (i.e., calculated using equations 24 and 26) whose significance is calculated to 98.97%. One complication with the AERONET data is that, apart from a significant signal, it is difficult to quantify exactly how changes in the Angstrom exponent translate into changes in various CCN characteristics. 5.2 Forbush Decreases and SSM/I Data The SSM/I measures the liquid water content over the oceans, and it is provided in a 1° × 1° grid. Each daily map is an average of the two adjacent days such that when the liquid water content is averaged over the oceans to give a time series, it is effectively a 3 day running average. Figure 5 shows the measured time series (equation 23) with three different statistical weights applied. In Figure 5a, the weights are based on FD strength (equation 17). The dotted lines on the left‐hand panel represent the significance calculated from the bootstrap samples displayed in equation 25. The corresponding integrated signal as defined by equation 24 from day t 1 =3 to t 2 =13 is shown on the right‐hand panel (red line). A different time interval than for the AERONET data is chosen since it is expected that the response of the clouds occurs later than that of the aerosols (see section 7.1). Each point (black dot) corresponds to one integrated response out of the 104 bootstrap samples. The significance of the integrated response is 99.68%. Figure 5b is similar to Figure 5a but assumes weights from equation 18. Here the integrated response is significant at the 99.99% level. Figure 5c is similar to Figure 5b except that it considers the five strongest FDs after 1998 (equation 18), thus allowing a comparison with the AERONET data. The integrated response here is significant at the 99.92% level. Finally, Figure 5d corresponds to all 26 FDs (equation 19) with equal statistical weights. The integrated response is here significant at the 95.72% level. For all the weights and time intervals considered, a significant reduction in SSM/I liquid water content of clouds over the oceans is seen in the days following FDs. It is also seen that the significance decreases as weaker FDs are added. This can be expected since weak FDs add significant noise without a significant signal. Figure 5 Open in figure viewer PowerPoint SSM/I data of liquid water content over the oceans. (left column) Superposition of SSM/I data using the strength of the FDs as (a) statistical weights (equation 17 ) and (b) equal weights for the five strongest FDs (equation 18 ). (c) Similar to Figure 5 b, but using only the five strongest FDs after 1998, and (d) uses equal weights for all FDs (equation 19 ). The red curve is a 3 day smoothing of the daily data. The hatched area denotes the temporal interval used to integrate the response (3–13 days). The dotted lines are 1, 2, and 3 standard deviations obtained using the MC bootstrap analysis. (right column) The bootstrap samples of the integrated signal. The abscissa is the bootstrap sample number and the ordinate is the size of the integrated bootstrap signal using the different statistical weights corresponding to the different rows. The red lines mark the size of the real FD signal with the corresponding weight. The achieved significance levels from top to bottom are 99.68%, 99.99%, 99.92%, and 95.72%. 5.3 Forbush Decreases and ISCCP Data The data from the Internationale Satellite Cloud Climate Programme (ISCCP) [Rossow and Schiffer, 1991] covers the period 1983–2006. Here the D1 data are used, which have global coverage and a 3 h temporal resolution (averaged into daily data for this study). Of the parameters offered by the ISCCP, we use both the IR‐detected total clouds and the IR‐detected low, middle, and high clouds. The IR detection of clouds is used since it has the lowest intrinsic noise compared to the visual channels; moreover, the observation of clouds over the oceans is more accurate and therefore used in the following [Brest et al., 1997]. Figure 6 shows the total IR‐detected cloud fraction over the oceans for the three statistical weights. Figure 6a shows the response where the individual FD time series have been weighted with the strength of the FDs. Figure 6b shows the five strongest FDs, while Figure 6c covers the five strongest FDs after 1998, which allows the comparison to the similar figures of the AERONET data and MODIS. Finally, Figure 6d assumes even statistical weights for all 26 FDs. In each of the panels the achieved significance levels for the integrated signal (3–13 days) is stated in the upper left corner. Figure 6 Open in figure viewer PowerPoint Statistics of the total IR‐detected ISCCP clouds over the oceans. (a) The variation in the total IR‐detected clouds averaged over the oceans using the strength of the FDs as statistical weights. (b) The five strongest FD (see, Table 1 ), (c) the five strongest FDs after 1998 (see, Table 1 ). (d) Based on even statistical weights for all the 26 FDs. The black curves denote the change in the response of the total IR cloud fraction, while the red curve is a 3 day smoothed version of the black curve. The hatched areas denote the temporal interval used to integrate the response (3–13 days), with the achieved significance levels stated in each panel. The dotted lines are 1, 2, and 3 standard deviations obtained using the MC bootstrap analysis. Figure 7 depicts the signal in all, high, middle, and low IR ISCCP clouds. It is seen that the signal is strongest for the total cloud product, whereas the signal is weaker in the individual cloud layer data sets. The reason for this could be that ISCCP satellites only detects the top of the clouds, and variations in, for example, lower clouds are disturbed by variations in the cloud layer above. This problem will be further discussed in the following section where we analyze the MODIS data. Figure 7 Open in figure viewer PowerPoint Statistics of IR‐detected ISCCP cloud fractions averaged over the oceans for the five strongest FD after 1998. The red curve is a 3 day smoothing of the daily data. (top left) All IR‐detected clouds. (top right) High clouds. (bottom left) Middle clouds. (bottom right) Low clouds. The hatched area denotes the temporal interval used to integrate the response (3–13 days) and the achieved significance levels are stated in each panel. The dotted lines are 1, 2, and 3 standard deviations obtained using the MC bootstrap analysis. Note that there is a discrepancy between Figure 7a and Figure 1d from Svensmark et al. [2009] which also depicts total IR‐detected ISCCP clouds. This is because Figure 1d from Svensmark et al. [2009] only includes data between 40° northern and 40° southern latitude. 5.4 Forbush Decreases and MODIS Data From “MOD08_D3” (the data product can be found on the MODIS website) we chose the parameters listed in section 2.3. The chosen parameters and their MODIS names are shown in Table 2. The daily average maps given by the data product cover all of Earth, and from these maps we produce global daily means for each of the parameters. The CCN product is ocean only due to difficulties with retrieval of aerosol counts over land [Levy et al., 2010]. Details on how the CCN product is derived can be seen in Remer et al. [2005]. Table 2. Parameters and Their Name in the “MOD08_D3” Data Product Parameter MODIS “MOD08_D3” Parameter Name ε (‐) Cloud_Effective_Emissivity_Mean CCN (108 cm−2) Cloud_Condensation_Nuclei_Ocean_Mean τ (‐) Cloud_Optical_Thickness_Liquid_Mean LWP (gm−2) Cloud_Water_Path_Liquid_Mean CF (‐) Cloud_Fraction_Liquid_Mean R eff (μm) Cloud_Effective_Radius_Liquid_Mean Figure 8 displays the result of the statistical test in the case of the five strongest FDs, in the period after MODIS data are available, starting with the event on 16 July 2000 (note that there are no strong FDs between 1998 and 2000, so these are the same FDs as in the AERONET analysis). The panels show the time series as defined in equation 16 for each of the above parameters using the weights in equation 18. The decreases in ε, τ, LWP, and CF as well as the increase in R eff extend beyond the 95% level. Figure 8 and Figure 1c of Svensmark et al. [2009] both depict the liquid CF as measured by MODIS but differ slightly in shape because Svensmark et al. [2009] treated “CF =0” values as missing data; an issue that has been corrected in the present work. Figure 8 Open in figure viewer PowerPoint ε, τ, CF, CCN, LWP, and R eff averaged for the strongest five Forbush decreases from year 2000 and onward (Table σ significance levels. The hatched area is the temporal range used to integrate the response, and the achieved significance levels (ASL) for this measure is stated on each figure (see section MODIS global daily means of the parameters, CF, CCN, LWP, andaveraged for the strongest five Forbush decreases from year 2000 and onward (Table 1 ). The black curve is the response in the cloud parameter, while the red curve denotes a 3 day smoothed version of the black curve. The dotted lines are 1, 2, and 3significance levels. The hatched area is the temporal range used to integrate the response, and the achieved significance levels (ASL) for this measure is stated on each figure (see section 4 ). Six days after the FD minimum, the parameters LWP, τ, and CF reach their minimum values, compared to ε which takes 8 days to reach minimum value. At day 9 the minimum for CCN is reached; however, it is less than 80% significant. The effective radius (R eff ) reaches its maximum after 11 days. The days when the parameters reach their extremum values when the average of the five strongest individual events is used, are shown in column 5 of Table 3. The ASL of integrated responses, using equation 28 from days 3 to 13 following the FD, is displayed in each of the figures for the six parameters in Figure 8. It is seen that all FD signals are significant at the 98% level or better, except for the CCN parameter, which will be discussed in section 6. Table 3. List of Quantities for the Six MODIS Parameters Using Data From FDs 1–5 Parameter Reference Level±σ ΔP meas (%) ΔP der (%) Extremum (days) ε (‐) 0.686±0.003 −1.29 −1.7±0.78 [30] 7.7±4.5 CCN (108 cm−2) 2.60±0.06 −3.32 −2.5±5.3 [32] 6.1±4.0 τ (‐) 11.09±0.12 −2.87 −3.7±1.5 [31] 8.1±4.5 LWP (gm−2) 108.60±1.11 −3.05 −2.2±1.6 [31] 8.5±4.1 CF (‐) 0.277±0.004 −5.53 – 9.5±4.3 R eff [μm] 16.95 ± 0.07 0.71 −0.19 ±2.1[31] 6.9 ± 4.3 We quantify the connection between Forbush decrease strength and parameter response for each of the MODIS parameters, in the same way as was done for the ISCCP (low clouds), MODIS (liquid clouds), AERONET (Angstrom exponent 330–440 nm), and SSM/I (liquid water path) data sets in Svensmark et al. [2009]. As the independent variable, we use the percentage change in ionization as listed in Table 1 for all of the FD events. The corresponding dependent variable is defined as a percentage change (ΔP w(i) ) in P w (t) i , where i refers to one of the six cloud parameters, and ΔP w(i) is the difference between the mean and the extremum value between days 0 and 15 of the FD. Minima in the data sets were used for all parameters except for the effective radius (R eff ) where maxima were used. Figure 9 shows the resulting scatterplots for all of the six parameters, together with least squares linear fits and LOESS regression. FD 1, farthest to the right on the plots, is an outlier which can perhaps be understood since the FD happened at the same time as the Halloween event of 2003, which contained eleven X class solar flares during the 18 days of the event [Woods et al., 2004]. Figure 9 Open in figure viewer PowerPoint t test, this does not mean that the relationship is necessarily linear. To illustrate this we have also added a LOESS regression (red lines) using a second‐order polynomial and a 0.8 fraction of the total points in each local regression. Also note that equivalent figures for the ISCCP, SSM/I, and AERONET data have already been published in Svensmark et al. [ 2009 The Forbush decrease strength (from Table 6 ) and response in six MODIS cloud parameters. The black lines are weighted linear trends of the data points. Slope values and standard deviations of the slope are written on each plot. The broken lines show the same trend, except that the data point for the extreme Halloween Event (FD 1) has been excluded. Although linear fits are used to permit atest, this does not mean that the relationship is necessarily linear. To illustrate this we have also added a LOESS regression (red lines) using a second‐order polynomial and a 0.8 fraction of the total points in each local regression. Also note that equivalent figures for the ISCCP, SSM/I, and AERONET data have already been published in] using the same list of FDs. Student's t test is used in order to determine if the fitted slopes are statistically different from zero. The significance of the slope of ε is 98% (99% excluding FD 1), for LWP it is 95%, for CF 95% (90% excluding FD 1), R eff 95% excluding FD 1 and insignificant when including it, and for τ it is 95% (99% excluding FD 1). Regarding CCN the slope is significant at the 90% level when FD 1 is excluded, and it is insignificant when including it. If the reverse analysis is carried out, i.e., looking for minima in R eff and maxima for the other parameters, we find that the probability for the slopes to be different from zero to be insignificant, which indicates that the result is not due to a symmetric increase in the level of fluctuations. So when the GCR influx decreases, it is followed by a decrease in ε, CF, LWP, and τ and possibly an increase in R eff . Note that this result does not necessarily imply that the relation is strictly linear, it merely suggests that there is a connection between the strength of the FD and the cloud parameter response. The LOESS results also shown in Figure 9 illustrate that the relation is not completely linear. The results in the last part of section 5.3 on the ISCCP data do not provide a clear answer to where in the atmosphere the effect can be seen. Returning to this question with the MODIS data, we now look at the division between ice clouds fraction and liquid cloud fraction. Figure 10 demonstrates that the signal is almost exclusively seen in liquid clouds. MODIS uses a number of near‐infrared and visible channels to determine the cloud thermodynamic phase, i.e., if a cloud region consists of liquid or ice particles [Platnick et al., 2003]. The combined signal (ice + liquid + undetermined clouds) in Figure 10a shows a clear response to FDs. However, if the plot is divided into the subtypes ice and liquid clouds (Figures 10b and 10c, respectively), then it is clear that the response originates from the liquid cloud fraction, which is consistent with the strong signal also seen in Figure 7 (top left) in the ISCCP All IR‐detected clouds. Note that in the study by Marsh and Svensmark [2000] the response in clouds to the 11 year solar cycle was clearly strongest for low clouds. In contrast, the present study can only confirm that the short‐term FD response is seen in liquid clouds, possibly due to the reasons mentioned in section 5.3. Figure 10 Open in figure viewer PowerPoint Response in MODIS‐derived cloud fractions by type, to the five strongest FDs from year 2000 and onward. (a) Cloud fraction combined. (b) Ice cloud fraction. (c) Liquid cloud fraction. The red curve is a 3 day smoothing of the daily data. The hatched area denotes the temporal interval used to integrate the response (3–13 days). The dotted lines are 1, 2, and 3 standard deviations obtained using the MC bootstrap analysis. 5.5 Intercorrelation in MODIS Data A possible issue with multiparameter analysis is intercorrelation, which can lead to overestimating the significance of simultaneous signals in the investigated parameters. Causes for intercorrelation could be if parameters are used to derive each other or if some of the same optical channels on the instruments are used to measure multiple parameters. For instance, liquid water path is derived from optical thickness and effective radius [King et al., 1997, page 65], although effective radius and optical thickness are measured using different wave bands, as seen in Table 1 [King et al., 1997]. The intercorrelations between the six investigated MODIS parameters are shown in Figure 11. The highest correlation found is between liquid water path and optical thickness (r = 0.88) which could be expected since one is derived from the other, as described above. Looking at the panels of those two parameters in Figure 8, they also appear quite similar. The only other two parameters that have a correlation coefficient (r) above 0.5 is emissivity and cloud fraction. Figure 11 Open in figure viewer PowerPoint Intercorrelation of the six MODIS parameters. Below the diagonal from top left to bottom right, correlation plots are shown. Above the diagonal the correlation coefficients corresponding to the plots are shown. Data points above five sigma have been filtered from the data in order to avoid distortion from unphysical events (instrumentation errors). Furthermore, seasonal variations were removed by filtering the data with a Fourier filter, removing variations greater than 90 days. 5.6 Principal Components Analysis (29) Miller and Miller, 2000 To combine the information from several data sets we perform a principal components analysis (PCA) on the eight cloud parameter time series, namely, the six MODIS parameters, the ISCCP total IR clouds over oceans, and SSM/I. The PCA procedure constructs an orthogonal transformation matrix, defining (up to) eight new time series as linear combinations of the original eight as projections onto the new basis. The first principal axis is chosen to account for the highest possible amount of variance in the system, and the projection of the time series data along this direction referred to as the first principal component (PC). (The second principal axis accounts for the second most amount of variance, projections onto it provides the second principal component, and so on). The first PC can be interpreted as a measure of the total change in the cloud system during the five strongest events in the period starting with the FD on 16 July 2000. To use the PCA procedure, all parameters were normalized by subtracting the mean value, removing the linear trend, and finally dividing by the standard deviation of the time interval. The first PC iswhere the numerical coefficients are the first eigenvector of the correlation matrix which can be seen as the amount of variance in the time series with multiple variables [, chapter 8]. The six MODIS parameters can be seen in Figure 8 , the ISCCP parameter in Figure 7 (top left), and finally the SSM/I parameter in Figure 5 c. Figure 12 presents the first principal component over the 36 day period. The variance is determined by Monte Carlo sampling using random time series for each of the eight parameters and then calculating their first PCs follow. The standard deviation is calculated from these realizations. A clear minimum is seen at days 5 and 6, which demonstrates that a signal in the cloud parametric system is simultaneously found in all (or most) parameters. Integrating the response from days 3 to 13 and comparing with Monte Carlo realizations, it is found that none of the 104 realizations can display a similar large integrated signal, and the ASL > 99.99%. This is the case even if the absolute values of the area is used. It is also important to note that the sign of the components of the first PC is found to have the R eff parameter with a sign opposite to all the other parameters. This is precisely predicted by cloud physics (see next section), and the probability of having this particular sign relation for the eight parameters from a random realization is even more unlikely than any of the isolated calculated significance levels. The CCN parameter from the MODIS data contributes with low variance to the first PC and its sign is wrong; however, as will be discussed in the next section, the CCN signal is expected to be unobservable due to the high noise level of the MODIS CCN data product. Due to the intercorrelation between liquid water path and optical thickness in the MODIS data, it could be argued that only one of these parameters should be used, but we include both for the sake of completeness, and since they are also both included in the MC analysis the significance should not be affected artificially by this. In summary, the PCA appears to strengthen the conclusion that the cloud system is actually disturbed when FD events occur. Figure 12 Open in figure viewer PowerPoint The first principal component (based on the six MODIS parameters, ISCCP total IR clouds over oceans and SSM/I) averaged over the five strongest events in the period starting with the FD on 16 July 2000. The dashed lines show 1 and 2 standard deviations. The hatched area from days 3 to 13 is the interval over which first principal component is integrated. The achieved significance level from the integrated signal here is >99.99%, meaning that none of the 104 Monte Carlo realizations gave an absolute result of the same size or larger.

6 Estimated and Observed Changes in Cloud Parameters We now turn to the physics of the measured changes in the observed parameters and calculate whether the changes are within the expected range. Based on observations of the energy entering the oceans over a solar cycle, a peak to peak variation of 1–1.5 W/m2 is found in the radiative budget [Shaviv, 2008; Howard et al., 2015]. This change has been associated with an absolute change of ≈2% in low cloud fraction (corresponding to an ≈5% relative change in low cloud fraction). Since the effect of the largest FDs is only slightly smaller than a solar cycle variation, as seen in Figure 3, the expected variation during the strong FDs is also of the order of 1–2 percentage points (or possibly smaller, since the effect in total liquid clouds may be smaller than that in low clouds alone). The question is how large such a signal is compared with the intrinsic noise of the various data sets. Stephens[ 1978 Hobbs [ 1993 (30) Column 2 of Table 3 lists the reference levels of the six MODIS parameters, and the mean percentage changes during the five largest events are shown in Column 3. We define the reference level according to equation 4 , and the percentage change relative to the reference using the extremum value within days 3 and 13. Now we want to check if the magnitude of changes in the MODIS parameters are internally consistent and if we would expect to see a signal at all, considering the noise in the data. To test this we used a series of equations from] (equation 30 ) and Chapter 2 of] (equations 31 ): (31) (32) Here a 0 is a scaling parameter which is found from the MODIS data by using reference levels for ε and LWP and then solving for a 0 , ρ is the density of water (1000 kg m−3), and N c is the droplet column density where CCN is used as an approximation. Using these equations we can then take a given change in one of the MODIS parameters and calculate the magnitude of change we would expect in some of the other parameters, the derived percentage change (ΔP der ). P der is found for each parameter using the reference levels (as seen in Table τ, is found for each parameter using the reference levels (as seen in Table 3 ) of the parameters they can be calculated from along with corresponding changes during a FD. For example, for The calculated ΔP der for all parameters is shown in Table 3. Note that the measured percentage drop of −2.87% for τ lies within the uncertainty of the derived value. Using equation 32 the derived change for N c is found to be ΔN c,der =−2.5 ± 5.3%. With the assumption that CCN changes as N c , ΔCCN der is within one sigma of the CCN data from Figure 8. The lack of significant CCN signals in Figures 8 and 9 can possibly be explained by the fact that the expected change is smaller than the noise. Similarly, ΔR eff, der is also contained in the noise. For LWP, ε, and τ the observations and derived parametric changes appear to be consistent, and this remains the case for parameters where several of the equations ((30)–(32))-((30)–(32)) could be used to derive ΔP der . Table 3 summarizes the results. Turning to the observed extremum variation in the other data (SSM/I, ISCCP, and AERONET), Table 4 summarizes the observed mean values and the extremum change in percent following the five strongest FD. Also shown are the standard deviation of the fluctuations in percent of the mean value. The signal‐to‐noise ratio is estimated from these values, from which one can see that the SSM/I has the largest signal‐to‐noise ratio, whereas the ISCCP IR low clouds have the smallest ratio. Table 4. Mean, Noise, and Signal Levels for the Four Data Sets Using the Five Strongest FDs Evenly Weighted Time Seriesa Observational Platform Parameter Mean Value Noise Level Noise Level (%) Signal Size (%) Signal/Noise MODIS Liquid cloud fraction 0.355 (‐) 0.005 (‐) 1.5 3 2.0 ISCCP All IR cloud fraction over oceans 0.64 (‐) 0.006 (‐) 0.9 2.3 2.6 ISCCP Low cloud fraction over oceans 0.32 (‐) 0.008 (‐) 2.5 2.5 1.0 SSM/I Liquid water path 0.090 (kg/m2) 0.001 (kg/m2) 1.1 3.3 3.0 AERONET Angstrom exponent 1.25 (‐) 0.05 (‐) 5 8 1.6

7 Discussion The fundamental question addressed here and in previous analyses is whether a real physical response exists in the aerosols and cloud parameters in the days following Forbush decreases. By ranking FD events according to their impact on the ionization in the lower atmosphere (Table 1), it was possible to test not just the atmospheric response to FD but also the size of the response as a function of the FD strength. The test for significance was performed using a Monte Carlo bootstrap statistical procedure defined in section 4. A significant signal was found in all of the four independent data sets in the days following the FD minimum (AERONET, SSM/I, ISCCP total IR clouds, and all MODIS parameters except for the MODIS CCN parameter, which, however, is expected based on its inherent noise level). Second, applying the Monte Carlo bootstrap statistical method by integrating the response in the days following the FDs leads to high statistical significance of the observed responses. Achieved significance for the five strongest FDs after 1998 are, AERONET: 98.97 (Angstrom exponent), SSM/I: 99.92% (liquid water content), ISCCP: 99.98% (All IR cloud fraction), and MODIS CF: 99.90% (liquid cloud fraction). The integrated response automatically addresses autocorrelations in the data sets. Equally important is that the numerical changes in the different cloud parameters are found to be consistent with the expected changes in the physics of aerosols and clouds as discussed in section 6. The consistent chain of reactions is fewer cosmic rays less ionization less aerosol nucleation, fewer formed CCN fewer cloud droplets larger cloud droplets [Boucher et al., 2013], decrease in cloud fraction, cloud optical thickness, and in cloud emissivity. Finally, since the droplets are larger, their rainout is more effective and it is consistent with the reduction in liquid water content. Furthermore, the atmospheric responses are found to scale with the strength of the FD, as seen in Figure 9. One sees that the signals are weaker for the less powerful FDs, suggesting that the response of the weak FDs is dominated by noise as could be expected. A likely reason is that the ionization and other average cloud parameters are fluctuating due to meteorological changes, and thus mask the effect of small FDs. The errors in the estimation of FD strength can also play a part in masking an effect from the weaker FDs. 7.1 Signal Delay The time between the FD and the extremum of the response in the cloud parameters is referred to as the signal delay. For the cloud data sets investigated in this study the signal delay is consistently between 6 and 11 days (see also Svensmark et al. [2009]). If the mechanism behind the changes is due to effects of ionization on the nucleation and growth of aerosols, it is expected that there would be a substantial signal delay as several processes would have to happen before a signal would be visible in the cloud data. The first step is that freshly nucleated aerosols (about 1 nm in size) would have to grow to CCN sizes, which is up to ∼100 nm depending on atmospheric conditions. A growth rate for aerosols in the order of 1 nm h−1 is not uncommon [Kulmala et al., 2004], which can explain a signal delay of several days, such as is seen in the AERONET data (section 5.1). The second step is that the CCN‐sized aerosols have to activate in order to grow into cloud drops, and furthermore. the new drop concentration (and size distribution) will affect the clouds through, e.g., rain out if there is fewer (and therefore larger) cloud drops. Adding up these steps, it is thus not unreasonable to observe a signal delay of about a week. We thus choose days 3 to 13 as the time period where a signal can be expected, for the cloud data sets. Varying this by 1 to 2 days in either direction can cause some change in the signal, but not in a way that changes the overall conclusions. Note that if the above response is caused by a decrease in aerosol nucleation and subsequent response in CCNs, it would seem to go against the common theory as implemented in large global aerosol models [e.g., Snow‐Kropla et al., 2011]. In such models, any small change in the formation of new aerosol particles (caused by, for example, a change in ions) vanishes before growth to CCN sizes can be reached. Thus, if ions are responsible for the observed changes in the present analysis, then ion effects in current atmospheric models are underrepresented. 7.2 Ultra Violet and Total Solar Irradiance Although the above results are consistent with a cosmic ray/atmospheric ionization link, effects on cloud cover could also arise from sensitivity to changes in the UV and/or the Total Solar Irradiance (TSI). Traditionally, it has been difficult to attribute cloud changes to one process or the other [Boucher et al., 2013; Laken et al., 2011], so in this section we investigate a possible UV and/or TSI effect. 7.2.1 Ultra Violet Irradiance (33) hν represents a photon, initiates most of the sulfuric acid production since O(1D) produces the OH radical, which reacts with SO 2 to form sulfuric acid. The photolysis constant can be calculated from (34) λ is the wavelength, t time, S(λ,t) is the solar spectrum at the top of the atmosphere, τ(λ) is the transmission through the atmosphere, σ(λ) is the cross section, and finally, Y(λ) is the quantum yield function. Figure DeLand and Cebula, 2008 for the three strongest FDs since 2003. Notice that the relative variation is small, in particular, when compared to the typical variance shown as the dashed lines. Changes in the UV range are important for photochemistry in the atmosphere, and they could therefore influence aerosol formation and cloud microphysics. As a consequence, it is not unreasonable that a simultaneous change in relevant UV wavelengths during FDs could explain the change in aerosols and clouds. This is especially so as one of the most important trace gases for aerosol nucleation and growth is sulfuric acid. Thus, a decrease in the UV frequencies concurrent with the FDs would lead to less sulfuric acid, fewer nucleated particles, and slower growth. The following photochemical reactionwhererepresents a photon, initiates most of the sulfuric acid production since O(1D) produces the OH radical, which reacts with SOto form sulfuric acid. The photolysis constant can be calculated fromwhereis the wavelength,time,) is the solar spectrum at the top of the atmosphere,) is the transmission through the atmosphere,) is the cross section, and finally,) is the quantum yield function. Figure 13 (top) demonstrates that the range of wavelengths where the photoreaction occurs in the lower part of the atmosphere is quite narrow. Using observational data of the solar spectrum from either composite solar spectral irradiance covering the time period 8 November 1978 to 1 August 2005 [] in the wavelength range 120–400 nm or the Solar Radiation and Climate Experiment (SORCE) covering the temporal period 14 April 2003 to 24 August 2015 ( http://lasp.colorado.edu/home/sorce/data/ ), it is possible to study the changes in the above photochemical reaction during FDs. The results all show that the variations are of the order of a few parts per thousand. Namely, they are too small to be important. As an example, Figure 13 (bottom) depicts the variation infor the three strongest FDs since 2003. Notice that the relative variation is small, in particular, when compared to the typical variance shown as the dashed lines. Figure 13 Open in figure viewer PowerPoint can occur in the lowest part of the atmosphere. The relatively narrow shape is a result of the atmospheric transmission τ(λ) cutoff at shorter wavelengths and the cross‐section σ(λ) cutoff at longer wavelengths. (bottom) Composite change in the photolysis constant for the reaction in equation J, the change is only of the order 2–3 × 10−3. The dashed lines are the 1σ variance. (top) The wavelength range where the reactioncan occur in the lowest part of the atmosphere. The relatively narrow shape is a result of the atmospheric transmission) cutoff at shorter wavelengths and the cross‐section) cutoff at longer wavelengths. (bottom) Composite change in the photolysis constant for the reaction in equation 33 , for the three largest FDs in the interval 2003–2006. Although there is a decrease in, the change is only of the order 2–3 × 10. The dashed lines are the 1variance. Moreover, the decrease in UV seen in Figure 13 (bottom) is not seen at wavelengths shorter than approximately λ = 280 nm, but an increase. Figure 14 displays the superposed change in UV at the wavelength λ = 250 nm centered on the three strongest FDs after 2003. Wavelengths shorter than λ = 280 nm do not penetrate into the troposphere, but UV with a wavelength λ = 250 nm has the largest absorption by ozone in the stratosphere. These shorter wavelengths are therefore important at the stratosphere. Figure 14 Open in figure viewer PowerPoint λ = 250 nm superposed and centered at the three strongest FD after 2003 (data from FD 2 and FD 5 from Table nσ variance, n = 1,2… . UV variation at= 250 nm superposed and centered at the three strongest FD after 2003 (data from FD 2 and FD 5 from Table 1 are not available in the SORCE data set). Notice that the relative change is of the order 0.8%. The maximum seen in UV is due to the rotation of active regions with a period of 27 days (see also Figure 15 ). Typically, an active region at maximum UV and a FD are closely correlated since a coronal mass ejection has to hit the Earth to generate a significant FD. The dashed lines are thevariance,= 1,2… . A clear maximum in UV is seen at days 1 and 2, minima at days −5 and 14, respectively. The origin of these changes in UV are caused by the approximately 27 day solar rotation of active regions on the Sun. That this is the case can be seen from Figure 15, which depicts the change in UV at the wavelength λ = 240.92 nm over the 400 day period, 14 April 2003 to 15 October 2015. Here the changes in UV seem to follow the approximately 26 day solar rotation of active regions on the Sun. Since the FDs have their origin from coronal mass ejections from active regions, the intensity of the UV at this wavelength is close to maximum as the active region passes the solar disk, as is exemplified by the strongest FD in Oct. 2003 (marked with the red diamond symbol). This FD did not impact the UV in any notable way, which can be seen from the subsequent three solar rotations where a maxima in UV are seen but without any FD superimposed. Figure 15 Open in figure viewer PowerPoint UV variation at λ = 245 nm as a function of time after 14 April 2003. The strongest FD event in October 2003 is marked with a red diamond symbol. Notice that the typical variation in UV is about 26–27 days and that the strong FD is not seen in the UV data. The blue symbols are maxima of 27 days rotation periods without any recording of a strong FD. The dates marked with the blue stars are 26 April 2003, 26 May 2003, 27 November 2003, 23 December 2003, and 21 January 2004. The question is now if such 26−27 day variations of typically 1 to 2% in UV in the stratosphere has an effect on clouds in the troposphere. In Figure 15 a number of 27 days UV variations are seen with variations equal to or larger than the 31 October 2003 event, whose maxima are marked with blue star symbols. The dates are 26 April 2003, 26 May 2003, 27 November 2003, 23 December 2003, and 21 January 2004. To answer the aforementioned question, Figure 16 displays the cloud fraction using “ISCCP all UV detected” data, where five time series are superposed over 36 days with day 0 being the date of maximum UV (shown as with the blue symbols in Figure 15). As can be seen from Figure 16, no statistical significant response is found (ASL = 83 %). It should be mentioned that since FDs are related to active regions, there is also a correlation between these active regions and cosmic ray variations. However, the dates used above for the UV maxima are not associated with any significant or strong FDs. Figure 16 Open in figure viewer PowerPoint n − σ variance, with integer n. Cloud fraction from the ISCCP all UV data superposed over the five dates shown in Figure 15 with maximum UV. The dashed lines denote thevariance, with integer The relative changes in tropospheric photochemical reaction rates simultaneous with strong FDs are of the order of parts per thousand. This should be compared with the typical simultaneous changes in clouds of ≈3%. A cause and effect from photochemistry is unlikely since it would imply an amplification factor of ≈30. UV wavelengths relevant to the stratosphere do show changes of the order 1–2% simultaneously with the strong FDs; however, similar changes in UV but on dates without strong FDs do not result in any significant cloud response. UV in the stratosphere is therefore unlikely to be the cause of the observed cloud changes during strong FDs. Summarizing the above, the data suggests the following: 7.2.2 Total Solar Irradiance Total solar irradiance (TSI) exhibits a decrease with a minimum occurring approximately 2 days before the FD minimum, as seen in Figure 17 for five strongest FDs. The observed decrease is of the order 1 W/m2, corresponding to a relative change of the order 0.0007—a relative change similar to the change in the photolysis constant. Such a change, which after distributing the energy over the Earth and taking the Earth's albedo into account, is of the order 0.2 W/m2, which is too small to make any significant impact on the atmosphere through any thermodynamic effects. It has been suggested that changes in the solar spectrum (mainly in the UV) can cause a warmer stratosphere that subsequently couples down to the troposphere over the 11 year solar cycle. But in the case of 0.2 W/m2 over a few days, due to the large heat capacity in the Earth's system, the change is so minute that any thermodynamic response is negligible. Second, a priori there is no mechanism that should relate the TSI with the observed changes in various CCN characteristics, such as their radius, since it depends on aerosol formation processes prior to the formation of clouds. Finally, similar to the argument used in the case of the photolysis constant, the relative variation in TSI is of the order 0.0007 and the corresponding changes in clouds are 0.03, which would imply an amplification factor of about 40 of the effect of TSI on clouds. This is unlikely. Figure 17 Open in figure viewer PowerPoint Variation of Total Solar Irradiance (TSI) around the mean for the five strongest FDs after 1998. The mean is 1365.6 W/m2 and the relative change is 0.001. The dashed lines denote the n − σ variance, with integer n. 7.3 The Global Electric Circuit An alternative theory about how Forbush decreases may affect the atmosphere involves the downward ionosphere‐Earth current density J z that is a part of the global electric circuit. For example, Tinsley and Deen [1991] reported impacts of FDs on winter storm vorticity. Another example is Kniveton et al. [2008], who examined changes in fair‐weather measurements of vertical electrical field at Vostok, Antarctica, and noted changes in ISCCP cloud data both at high latitudes and in the tropics. Although measurements of J z and the electric potential are highly dependent on regional and meteorological conditions, it is suggested from modeling that there is a systematic zonal change in J z following FDs and that such changes can influence cloud microphysics [Tinsley et al., 2007]. These processes may then influence aerosol concentrations several days later. J z has also been suggested to affect ice formation [Tinsley et al., 2007]. In the present work we do not see any effect in ice clouds but cannot rule out a J z effect in the liquid phase. 7.4 Comparison With Previous Studies We can now use the results in this study to explore why some previous studies did not find significant responses to FDs. Calogovic et al. [2010] examined six FDs from the ISCCP IR low clouds data set and concluded that there is no effect to be found. Based on the present study, we note the following likely reasons as to why this conclusion was reached: First, ISCCP IR low clouds have a fairly weak signal in contrast to ISCCP all IR clouds as can be seen by comparing Figure 7 (bottom right) with Figure 7 (top right) or the signal‐to‐noise values of Table 4. As discussed previously, this is likely due to the satellite view of low clouds being obscured by clouds at higher levels. Second, the selected FDs in Calogovic et al. [2010] all rank low in our list. Table 5 compares the FDs used by Calogovic et al. [2010] with the ranking in Table 1. These choices make it difficult to observe a signal. Using the same six FDs and data from SSM/I, which has the highest signal‐to‐noise ratio, and applying the procedure in section 4, we find a result with a marginal 95.7% significance. Table 5. List of the Six FDs Used in the Study of Calogovic et al. [ ] and Their Corresponding Ranking According to the List in This Studya Calogovic Et Al. Rank Order by Strength According to Table 1 1 13 2 8 3 5 4 10 5 24 6 15 Kristjánsson et al. [2008] used MODIS data to examine the means of 22 FDs. For comparison, only 13 FDs in the same period of 2000–2006 were used in the present study, and most of the 22 FDs resulted in minor changes in ionization, with the result that the mean signal was obscured by the meteorological noise. This can be the reason Kristjánsson et al. [2008] found that there was no signal. When they looked at the six strongest events on their list (these events are also high on the ranking list in this study) their signal improved. Sloan and Wolfendale [2008] used ISCCP data but focused in part on monthly averaged data. As shown in the present study (section 7.1) and in Svensmark et al.[2009], the extremum of the effect occurs about a week after the FD minimum and last only a few days; monthly averages are unlikely to show a signal as the signal will be obscured by noise. Laken et al. [2009] looked at MODIS data and considered a longer time series than 20 days (±40 days). They state that the maximum time for an aerosol particle to grow to CCN size is 2 days and conclude that a delay of about 5 days cannot be justified by any known process. However, a growth time of 5–7 days is in agreement with observations, as discussed in section 7.1. Moreover, they assess that the fluctuations observed during FDs are more likely noise than a signal related to cosmic rays. However, they did not carry out a statistical study; e.g., there are no error bars/confidence intervals on their figures. The study by Svensmark et al. [2009] was criticized for estimating the significance of FD signals based on an average variance over a 36 day interval, instead of a Monte Carlo estimate. This problem is avoided by the present approach.

8 Conclusion By ranking the strength of Forbush decreases according to their expected impact on the ionization of the lower atmosphere, a FD strength‐dependent response could be investigated. A Monte Carlo bootstrap statistical procedure was defined which either considered the time series or the integrated response over the days following the FD minimum and was based on the expected growth time of aerosols. These statistical tests allowed for a calculation of an achieved significance level and so could provide information on the likelihood that there was a response in atmospheric aerosol and cloud parameters following FDs. Four independent atmospheric data sets were used: (1) AERONET data using the Angstrom exponent in the wavelength range 330–440 nm. This parameter is related to the fine aerosol fraction. (2) SSM/I data measuring the cloud liquid water content over the oceans, (3) ISCCP data using IR detection of low, middle, high, and total cloud fraction, and finally, (4) MODIS data which allowed the study of a number of cloud microphysical parameters simultaneously. The parameters were cloud effective emissivity, cloud optical thickness, liquid water cloud fraction, column density of the cloud condensation nuclei, liquid water path, and finally, liquid cloud effective radius. Responses (>95%) to FDs are found in almost all parameters of the analyzed data sets: AERONET: Angstrom exponent, a measure of aerosols and cloud condensation nuclei [CCN] changes; SSM/I: liquid water content; ISCCP: total, high, and middle IR clouds above oceans; MODIS: cloud effective emissivity, cloud optical thickness, liquid water, cloud fraction, liquid water path, liquid cloud effective radius. ISCCP low IR clouds above oceans are only significant at a 93% level. In this connection it was observed that MODIS liquid cloud fraction also has a clear response in contrast to MODIS ice cloud fraction, which indicates that the effect is mainly in liquid clouds. Since the total UV‐detected clouds by ISCCP show a strong response, it suggests that this parameter is mainly influenced by liquid clouds. In contrast, variations in ISCCP low UV‐detected clouds are influenced by overlap of clouds at other heights, since the satellites have only seen the top layer of clouds and is the likely reason for the small significance of this parameter. One MODIS parameter, column density of the cloud condensation nuclei, is found to be insignificant which is expected based on its high noise level; however, independent data from AERONET do show a response. A positive, nonzero relation between the strength of the FDs and the size of the responses is found in all data sets. Changes in UV or in TSI were found to be unlikely to explain the observed responses in clouds or aerosols following FDs, since that would require an amplification factor of 30–40. Responses are found in the cloud microphysical parameters in the days following Forbush decreases. The sign and size of the response in all the parameters are consistent with changes derived from cloud microphysics. The size of the responses are of the order of a 2% change in cloud fraction for the strongest FDs. A correlation between the magnitude of the FD events and the effect on aerosol/cloud physics has been found in all data sets (AERONET, ISCCP, MODIS, and SSM/I). The signs of the responses are as expected from a cosmic ray effect on cloud microphysics. We therefore conclude the following: These results show with high confidence that there is a real impact of Forbush decreases on cloud microphysics. The suggested causal chain of reactions responsible for the observed correlations begins with a solar coronal mass ejection resulting in a FD with fewer cosmic rays less atmospheric ionization less aerosol nucleation fewer formed CCN fewer cloud droplets larger cloud droplets, decrease in cloud fraction, cloud optical thickness, and in cloud emissivity. Finally, since the droplets are larger, removal by rain is more likely and is consistent with the reduction in liquid water content. We note that a J z mechanism cannot be ruled out. In conclusion, the results support the suggestion that ions play a significant role in the life cycle of clouds.

Acknowledgments We are grateful to the late Nigel Calder for helpful comments in the early stages of this work. We thank Torsten Bondo for his initial work on part of the programs used in the data analysis. We thank professor of statistics Henrik Spliid for valuable discussions regarding the statistics. We thank the two reviewers for their thorough and constructive reviews which have helped improve the paper. We acknowledge receipt of the unpublished TSI data from the VIRGO Experiment on the cooperative ESA/NASA Mission SOHO (version d41_61_0803) from PMOD/WRC, Davos, Switzerland. We are also grateful to the Cosmic Ray Section, Solar‐Terrestrial Environment Laboratory, and Nagoya University who provided the muon data. Data can be obtained from http://www.stelab.nagoya-u.ac.jp/ste-www1/div3/muon/muon1.html We thank Tom Woods (CU LASP), Gary Rottman (CU LASP), and Giuliana de Toma (NCAR, HAO) for the SOLSTICE data and Mark Weber (U. Bremen, Germany) for the GOME data. MODIS data were obtained from http://modis.gsfc.nasa.gov. SSM/I data are produced by Remote Sensing Systems and sponsored by the NASA Earth Science REASoN DISCOVER Project. Data are available at http://www.remss.com. We thank the principal investigators and staff of AERONET (http://aeronet.gsfc.nasa.gov) for establishing and maintaining the sites used in this investigation. Martin B. Enghoff thanks the Carlsberg Foundation and the Danish Council for Independent Research/Natural Sciences for financial support.