Review of Sulo et al “Measurement Report: Increasing trend of atmospheric ion concentrations in the boreal forest”

The manuscript introduces an interesting result on increasing atmospheric ion concentration, especially as it has been shown that the aerosol particle concentrations are at the same time decreasing. The authors state that the manuscript gives the first proper assessment to atmospheric ions, which increases the value of this work if true. As other measurement sites are also starting to have long ion time series, this would be a nice reference point for comparing trends in different environments. The topic fits well to scope of ACP and I could recommend publishing it if the concerns below will be properly addressed.

Major comment:

My main concern is that the statistical methods are inadequately described and confusion in terminology raises concern that were the methods used properly and thus are the results valid. See my specific comments for details.

 Answer: We thank the reviewer for these important questions that will improve the clarity of our manuscript. We have answered to each specific comment below. We believe that the added discussion on the statistical methods and better explanations of the terminology have now clarified the methods and their uses in the paper and thus showed the validity of the results.

Specific comments:

RC: Page 2 lines 63-64: Equations are rarely seen in Introduction section, as it is meant to give more general background information. Consider restructuring the text.

Answer: We have moved the explanation of coagulation sink into Materials and Methods.

RC: Page 3, lines 89-90: Does the change of instrument location affect the radiation measurements? 

Answer: By visual inspection of the data there is no indication that it does. There does not appear to be a breakpoint in the ionization rate data for radon, the only radiation measurement available from the new site C as well as the old location. However, due to the sparsity of the data, this cannot be fully discounted either. We have added a mention of this to the text (Page 7: Line 165-168):

The sparsity of the in-situ ionising radiation measurements makes it difficult to detect breakpoints, but visual inspection revealed no breakpoints in the data. However, missing data periods in the radiation data can possibly hide breakpoints and changes due to the change in the measurement location cannot be conclusively dismissed. The condensation sink and meteorological data time series did not exhibit any visually noticeable breakpoints.

RC: Page 7, lines 161-164: It is invariably true that the presence of autocorrelation needs to be accounted for but prewhitening is not necessary the best method for that. PW loses information on the time series and thus sometimes it may cause a significant trend not to be detected (see e.g. Razavi&Vogel, 2018). If the trend can be assumed monotonic, then Sen’s slope is a proper choice and the significance for that can be tested without MK-test, and thus without prewhitening. This has been done e.g. in Leinonen et al. (2022) with bootstrap-based confidence intervals. Sen’s slope can be correct for autocorrelation with method by Kunsch (1989), if needed. If the trend cannot be assumed monotonic, then more proper methods should be used, like dynamic linear models (DLM, Laine 2020). With DLM, also seasonal variation and other cyclic structures or sources of measurement error can be accounted for.

Answer: We have applied the 3PW prewhitening method proposed by Collaud Cohen (2020) on those time series where the trend can be assumed monotonic. The 3PW method combines the traditional prewhitening method (PW) with the trend-free pre-whitening method created by Yue et al. (2002) and the variance-corrected prewhitening method presented in Wang et al. (2015). The Sen’s slope is then calculated on the variance-corrected trend-free prewhitened data, which leads to more accurate slope estimates than the other prewhitening methods. The MK-test used to test for significance of the trend is modified for seasonality as per Hirsch et al. (1982). We have now added a short description of the method into the text and more clearly point the reader to the Collaud Cohen (2020) article for more details. As per suggested by the reviewer, we have also added the dynamic regression trends calculated by a dynamic linear model as presented in Laine (2020). As can be seen from the modified figures 5 & 6 showing both these methods, the trends from the Sen’s slope and DLM are remarkably similar.

RC; Page 7, lines 165-166: How were the validity of assumptions for n-ANOVA tested?

Answer: The normality of the data was estimated visually from histograms. The data was deemed roughly normally distributed. Linear regression has been shown to be fairly robust even with violated normality assumptions (Knief and Forstmeier, 2018). All data was normalized to between 0 and 1 before applying multiple linear regression. We have added a Belsley collinearity test, which lead to the elimination of the cosmic ray ionization rate from the model. Additionally, upon review we noticed the test was utilizing an uncorrected version of condensation sink and we corrected this mistake as well in the revised manuscript. We rewrote chapter 3.2.3 to reflect the updated multiple linear regression results  and included descriptions of the assumptions and how they were met:



We defined our linear model as

N_{ion =} ωCS + βI_R + γI_g + δT + θRH + μI_CR + ε                                                                       (2)

where ω, β, γ, δ, θ and μ are model coefficients and ε is the error estimate.

Belsley collinearity test revealed that I_CR exhibited multicollinearity with RH. We therefore eliminated I_CR from the model and redefined it as

N_{ion =} ωCS + βI_R + γI_g + δT + θRH + ε                                                                       (3)

The model assumes that the inputs are normally distributed, independent and that the variance of the factors is roughly proportional.  To account for this, each factor was normalized by using the formula

x_norm = (x - min(x)/(max(x) - min(x))                                               (4)

where x_norm is the normalized data, x is the original data, min(x) is the minimum value of the data and max(x) is the maximum value of the data. The normalized data was then input into the model. The residuals of the models were normally distributed and we concluded that the model is usable for our analysis. The coefficients of the models are listed in Table A1.

Knief, U., Forstmeier, W. Violating the normality assumption may be the lesser of two evils. Behav Res 53, 2576–2590 (2021). https://doi.org/10.3758/s13428-021-01587-5

RC: Table 3: Write in caption what does the interval within brackets indicate (same thing with Figures 4 and 5)

Answer: We have added an explanation of what the values in the brackets mean:

“The values in brackets are the confidence intervals for the trend at a 90% level.”

RC: Section 3.2.3: Here is a slight confusion: Are you conducting here analysis of variance (ANOVA) or linear regression? From your results I would say the latter. Specify the model you are applying here and give the coefficients for the parameters at least in the appendix. Are the associations between variables assumed linear? Did you check how the residuals of your models look like? What is the time resolution of the data you are using in these analyses?

Answer: The reviewer is correct that our use of terminology is not clear. We estimate the amount of variance explained by each of our selected factors by fitting a simple linear multivariate regression model. We have specified the model type in the text and included the coefficients for the parameters in the appendix. The time resolution used is monthly medians.

RC: Figure 8: On what do you refer with R2? How were the individual R2 values derived? Are these sub-values from multivariable model (what is then the total R2?) or from separate bivariate models? If the model was multivariable, how was it constructed? Did you have all the predictors in the model at the same time, even if no significant? Did you check it for multicollinearity?

Answer: R2 here refers to the sum of squares for each variable divided by the total sum of squares of the multivariable model. We have added the total R2 of the models to figures 8,9 and A3-A6. All of the predictors were in the model at the same time, even if not significant. We have added a multicollinearity check and eliminated the cosmic ray radiation ionization rate from the model due to severe multicollinearity.

RC: Line 295: please elaborate what you mean about this. Do you mean seasonal variation or really variance (see definition for that)? how do you estimate that with monthly medians? How does your model for this look like? With properly selected models, you can take the seasonal variation and the long term trend account at the same time.

Answer: We have changed the wording to “seasonal variation” and further elaborated on our method in the text:

“We also analyzed the seasonal variation by using median concentrations calculated for each month. By calculating the median monthly concentrations of our data, we have a median yearly cycle for each factor, from which we can investigate how much the seasonal variation in ion concentration is explained by our model (Figure 10).”

References

Kunsch, H. R.: The Jackknife and the Bootstrap for General Stationary Observations, Ann. Stat., 17(3), 1217–1241, doi:10.1214/aos/1176347265, 1989.

Laine, M. : Introduction to Dynamic Linear Models for Time Series Analysis. in book: Geodetic Time Series Analysis in Earth Sciences, pp 139-156, doi: 10.1007/978-3-030-21718-1_4, 2020

Leinonen, V., Kokkola, H., Yli-Juuti, T., Mielonen, T., Kühn, T., Nieminen, T., Heikkinen, S., Miinalainen, T., Bergman, T., Carslaw, K., Decesari, S., Fiebig, M., Hussein, T., Kivekäs, N., Kulmala, M., Leskinen, A., Massling, A., Mihalopoulos, N., Mulcahy, J. P., Noe, S. M., van Noije, T., O'Connor, F. M., O'Dowd, C., Olivie, D., Pernov, J. B., Petäjä, T., Seland, Ø., Schulz, M., Scott, C. E., Skov, H., Swietlicki, E., Tuch, T., Wiedensohler, A., Virtanen, A., and Mikkonen, S.: Comparison of particle number size distribution trends in ground measurements and climate models, Atmos. Chem. Phys. Discuss. [preprint], https://doi.org/10.5194/acp-2022-225, in review, 2022.

Razavi S., R. Vogel: Prewhitening of hydroclimatic time series? Implications for inferred change and variability across time scales, Journal of Hydrology, Volume 557, 2018, Pages 109-115, https://doi.org/10.1016/j.jhydrol.2017.11.053.

Citation: https://doi.org/10.5194/acp-2022-392-RC1

The article by Sulo et al. describes and discusses the long term measurement of size resolved atmospheric ion concentrations in Hyytiälä/Finland for a time period of 16 years. The data set further includes the condensation sink, meteorological parameters and different parameters contributing to the ionization rate. Interestingly, a trend analysis shows that the ion concentrations increase over time. Atmospheric ions are relevant due to their ability to form new particles and grow them to larger sizes. Therefore, the findings discussed by Sulo et al. are important and fit well into the scope of the journal. I strongly suggest that the manuscript should be published in ACP after the points listed below have been considered.

RC: Recommendation for further analysis:

Given the unique data set, I think the authors should use it for a further analysis. They decided to show equation (1) but actually never use it. The equation could, however, be used to derive α (ion-ion recombination rate) from the measured n, CoagS and q by assuming steady-state conditions (dn/dt = 0). By plotting α as a function of time, significant deviations from the literature value (~ 1.7e-06 cm3 s-1) would indicate that a data point is not suitable for the further analysis (and should be filtered out). Or put differently, agreement between the literature value for α and the derived (calculated) values would indicate that the relevant data sets have been measured accurately.

Answer: We appreciate the reviewer’s suggestion but feel that a thorough investigation of the ion-ion recombination is out of scope for this manuscript, especially since it is submitted as “Measurement report”. More research is needed to properly quantify the ion-ion recombination rate from this dataset. The measurement accuracy of the ion concentrations are not the only reason why the recombination coefficient might differ from the literature values, especially since there is uncertainty also in the measured ionization rates and relevant loss processes. Franchin et al. (2015) showed that even in laboratory conditions, the experimentally determined recombination coefficient varied between ca. 1.5 - 9.7 × 10−6 cm3 s-1 depending e.g. on RH and temperature.

 However, we have done a quick calculation of the ion-ion recombination rate using our data set where all data for the ionization rates and coagulation sink are available. Because the ionization rates are theoretical maximum ionization rates and the uncertainties involved in them are comparatively large, our analysis should be not read as very conclusive. The median ion-ion recombination rate from our data was 4.8e-6 and the values varied within 2.2 x 10^-6 (5^th percentile) and 1.4 x 10^-5 (95^th percentile) cm³s^-1. This is somewhat higher than the value from Hoppel and Frick (1986), but interestingly, very close to the experimental values determined by Franchin et al. (2015). This indicates that the order-of-magnitude of the ion concentrations and ionization rates are correct, and as the reviewer suggests, is a good starting point for further analysis at a later stage.    

Franchin et al. Experimental investigation of ion–ion recombination under atmospheric conditions, Atmos. Chem. Phys., 15, 7203–7216, https://doi.org/10.5194/acp-15-7203-2015, 2015.

Minor points:

RC: Line 34/35: not a complete sentence, reformulate to, e. g. “…, typically the ions initially form from nitrogen and oxygen due to their high abundance in the atmosphere.”

Answer: We have edited the sentence to: “Atmospheric ions are produced via the ionization of air molecules (Rutherford, 1897) and primary ions are typically formed from nitrogen and oxygen due to their abundance in the atmosphere (Israël, 1970).”

RC: Line 40: please specify in the beginning that all diameters refer to “mobility diameters”

Answer: We have added a clarification to the text as follows: “...initially neutral particles of the same size (Tammet et al., 2013) or through ion-mediated nucleation (Hirsikko et al., 2011). All ion diameters mentioned here and further in the text are mobility diameters.”

RC: Line 43: “ions of this size are connected to …, snow fall or rain”; does this mean that ions initiate snow fall or rain or that they are produced from precipitation? Please clarify

Answer: We have clarified this by changing the sentence to: “Ions of this size are typically generated by atmospheric new particle formation (NPF), snow fall or rain (Manninen et al., 2010; Kerminen et al., 2018; Leino et al. 2016).  “

RC: Line 63: Please move equation (1) to section 2 and provide a value and reference for α.

Answer: We have moved this paragraph into section 2 and provided both a value (1.6e-6) and reference (Hoppel and Frick, 1986) for α.

Hoppel, W. A. and Frick, G. M.: Ion-aerosol attachment coefficients and the steady-state charge distribution on aerosols in a bipolar ion environment, Aerosol Sci. Technol., 5, 1–21, 1986.

RC: Line 98: “data are available”

Answer: We have corrected this to the correct form: “Figure 2: Time periods when data are available for...”

RC: Line 106: 0.82 nm are a rather small lower size limit, however, are smaller ions also existing, and if yes, how would their negligence in the measurement effect the outcome of the results?

Answer: The lowest channel of the BSMA (3.2-2.4 cm2/Vs ~ca. 0.82-1.07 nm) corresponds to roughly a mass range of ca. 30-100 Dal (Ehn et al. 2011), assuming a density of 1.66g/cm3. This means, that only the primary ions and some of the smallest molecular ions fall out of the lower size limit of the BSMA. Therefore, we do not believe this has a major effect on the cluster ion concentrations or the observed trends in this article, although for detailed studies on the ion balance or the transfer of charge from the primary ions to molecular ions and clusters, measurements of the primary ion size distribution and dynamics would be needed (Chen et al. 2016).

 We have commented on this in the discussion: “Also, for detailed studies on the ion balance and transfer of charge from primary ions to cluster ions, the measurements of primary ion size distribution below 0.8 nm would be needed (Chen et al. 2016).”

Chen, X., Kerminen, V.-M., Paatero, J., Paasonen, P., Manninen, H. E., Nieminen, T., Petäjä, T., and Kulmala, M.: How do air ions reflect variations in ionising radiation in the lower atmosphere in a boreal forest?, Atmos. Chem. Phys., 16, 14297–14315, https://doi.org/10.5194/acp-16-14297-2016, 2016.

Ehn et al. (2011). An Instrumental Comparison of Mobility and Mass Measurements of Atmospheric Small Ions, Aerosol Science and Technology, 45:4, 522-532, DOI: 10.1080/02786826.2010.547890

RC: Line 119/120: It is unclear what is meant by “if not necessarily its order of magnitude”.

Answer: We have removed the phrase to simply the text.

RC: Line 150/151: Do you have any ideas what the exact reason for the relation between noise and ageing is?

Answer: We hypothesize that the electrometer noise in the aspiration condenser increases as the instrument gets older and some particles are deposited into the condenser.

RC: Line 158: “particles”

Answer: We have corrected the word into plural.

RC: Line 163/164: Please provide references for the MK test and the Sen’s slope calculation. I also agree with referee 1 that the manuscript would benefit from further description of the analysis methods. 

Answer: We have included references to the tests used and added further descriptions of the analysis methods. Please note the answer to Reviewer 1 above for details.

RC: Line 180: It seem that the diameter range is not correct here as the mentioned ion concentration is too low for these small diameters.

Answer: The diameter range and the concentrations are taken directly from Manninen et al. (2009). The ion concentrations are much lower than the particle concentrations, which according to Manninen et al. (2009) are between 390 and 1290 cm-3. Note, that the size range 1.8-3 nm is mostly above the cluster ion band, and thus the ion concentrations in this size range are very low (outside periods of ion-induced nucleation) due to the low charging probability of aerosols at this size range.

RC: Line 195: “polarities”

Answer: We have corrected the term in the text.

RC: Table 3: Please replace the “-“ sign by the word “to” otherwise it can be confused with a minus sign.

Answer: We have replaced the n-dashes with the word “to” to make the table simpler to read.

RC: Figure 5: For the ionization rates: Does the sum of all three rates correspond to the value of q in equation (1)? What is the reason for the minimum in I_Radon between 2012 and 2014? I am not experienced in Radon measurements but I would assume that its value should not show a strong inter-annual variation. There also seem to be strong spikes for I_gamma, are these from known sources?

Answer: The sum of all three rates corresponds to the maximum theoretical ionization rate, but it does not necessarily correspond to the exact value of q in equation (1) because the ionization probability is not unity even with the requisite energy present. We cannot conclusively explain the radon minimum between 2012 and 2014. It is possible that instrument decay plays a role, as the instrument broke down in 2014 and was only fixed and returned to operation in 2017. The strong upward spikes in gamma were filtered out as they were the result of clear outliers in the original data, but the downward spikes are part of its annual variation clearly visible in the time series.

RC: Figure 6 and Figure 7: The figures are hard to read when printed out. It is also not clear to me what the red crosses indicate.

Answer: We have made the figures clearer and explained what the red crosses indicate. For example in Fig 7:



Figure 7: The annual variation of CS, T, RH, I_R, I_g, I_CR. The red line is the median concentration, and the blue box contains 50% (25^th to 75^th percentile) of all data points. The whiskers mark the location of the 95^th and 5^th percentile data points. The numbers on the x-axis are month numbers. The red crosses are outliers.

RC: Line 260: “about 50%”

Answer: This section has been rewritten as per the new statistical tests explained here and the slight change in the percentage of variability explained. We have added similar hedging as the reviewer suggests.

RC: Line 328/329: I do not think that this (BSMA/humidity effect) was mentioned before.

Answer: We have added a mention of this into section 2.2. (Data verificiation) and included a corresponding plot in the appendix (Fig. A2).

RC: Figure A1: What are the units on the y-axes in this plot? It seems that the smallest and largest size channels are missing (or at least the diameters do not agree with the size range mentioned in section 2.1.1).

Answer: We only use the size distribution up to 7 nm as we limit our analysis to the typical size range of small and intermediate ions (0.8 – 7 nm). We have corrected the size ranges to represent the actual size channels in the analysis and added labels to indicate the moving variance and its units (cm^-6) as well as size of the window used in calculating the moving variance.

Citation: https://doi.org/10.5194/acp-2022-392-RC2