Ensemble mean density and its connection to other microphysical properties of falling snow as observed in Southern Finland

In this study measurements collected during winters 2013/2014 and 2014/2015 at the University of Helsinki measurement station in Hyytiälä are used to investigate connections between ensemble mean snow density, particle fall velocity and parameters of the particle size distribution (PSD). The density of snow is derived from measurements of particle fall velocity and PSD, provided by a particle video imager, and weighing gauge measurements of precipitation rate. Validity of the retrieved density values is checked against snow depth measurements. A relation retrieved for the ensemble mean snow density and median volume diameter is in general agreement with previous studies, but it is observed to vary significantly from one winter to the other. From these observations, characteristic mass– dimensional relations of snow are retrieved. For snow rates more than 0.2 mm h−1, a correlation between the intercept parameter of normalized gamma PSD and median volume diameter was observed.


Introduction
Due to a variety of ice particle types and shapes, representation of winter precipitation in models (Woods et al., 2007;Morrison and Milbrandt, 2015) and in ground, airborne and satellite remote sensing retrievals (Sekhon and Srivastava, 1970;Matrosov, 1997;Wood et al., 2013) is a topic of continuous interest.Both models and retrieval algorithms rely on a prior knowledge of snowflake mass (or density), shape and fall velocity, which are typically expressed as functions of a characteristic particle size (Pruppacher and Klett, 1996).Furthermore, information on possible particle size distributions (PSDs) is also required.Even though some of the microphysical properties of ice particles are not independent, e.g., fall velocity can be computed from particle mass and shape (Böhm, 1989;Khvorostyanov and Curry, 2005;Mitchell and Heymsfield, 2005;Heymsfield and Westbrook, 2010), the remaining degrees of freedom are rather numerous.
Historically, measurements of snowflake properties have been carried out on a particle-by-particle basis (e.g., Magono and Nakamura, 1965;Locatelli and Hobbs, 1974;Mitchell, 1996).While we may still regard such measurements as the more precise and detailed, these studies are limited to a relatively small number of observed ice particles due to the Published by Copernicus Publications on behalf of the European Geosciences Union.
sheer amount of time needed for such experiments and corresponding data analysis.After the introduction of robust optical instruments capable of measuring particle size, shape and in some cases fall velocity, e.g., 2-D-video disdrometer (2-DVD; Hanesch, 1999;Schönhuber et al., 2007), particle size velocity (Parsivel) laser-optical disdrometer (Löffler-Mang and Joss, 2000;Löffler-Mang and Blahak, 2001), hydrometeor velocity size detector (HSVD; Barthazy et al., 2004), snow video imager (SVI; Newman et al., 2009) and multiangle snowflake camera (MASC; Garrett et al., 2012), continuous recording of ice particle properties became possible.It should be noted that, in comparison to surface-based observations, aircraft measurements have a much longer history in determining ice particle microphysical properties and were carried out in different types of clouds and climate regimes (Pruppacher and Klett, 1996).A typical limitation of automatic observations of ice particle properties, however, is that only a subset of needed parameters is directly measured.
By combining optical disdrometer observations with other measurements, e.g., by radar or precipitation gauge, physical properties such as mean snow density can be derived.Huang et al. (2010) have used a C-band weather radar observations of equivalent reflectivity factor, Z e , in combination with a 2-DVD to derive a snow density-dimensional relation and to infer more consistent Z e -snowfall rate relations.Another method for snow density retrieval is based on solving aerodynamic equations to derive particle mass from observed fall velocity and particle effective projected area as proposed by Böhm (1989) and applied by Hanesch (1999) and more recently by Szyrmer and Zawadzki (2010) and Huang et al. (2015).Brandes et al. (2007), hereafter referred to as B07, used a combination of a weighing gauge and a 2-DVD to derive a relation between mean bulk density and median volume diameter and to document relations between PSD parameters for Colorado winter storms.Their approach is similar to the one used by Heymsfield et al. (2004), who have combined aircraft PSD and ice water content observations to derive mean snow density and average mass-dimensional relations for ice particles.Albeit using slightly different definitions, both B07 and Heymsfield et al. (2004) derive effective ice densities for ensembles of ice particles, but there is a difference in terminology.Heymsfield et al. (2004) and many others have used the term (particle) bulk density to refer to the density of individual ice or snow particles defined as the ratio of mass of a particle with a size D to its assumed volume: ρ = ρ(D).In most of such cases, the word "bulk" is used to emphasize the inclusion of hollows within particles.The term "(mean) bulk density" is sometimes used also when referring to the mean density of an ensemble of particles representing the whole PSD, i.e., ρ = ρ(D 0 ) (e.g., B07), whereas Heymsfield et al. (2004) used the term "populationmean effective density".In this study we derive the volume flux weighted snow density, similar to, e.g., B07, and refer to it as ensemble mean density, ρ, to avoid possible confusion.This paper documents the connection between ensemble mean density and other microphysical properties of snow as observed in Southern Finland.Using the estimated ρ, average mass-dimensional relations characteristic to studied snowfall events are defined.In order to derive ensemble mean density, a method proposed by B07 was used.However, instead of a 2-DVD, a new generation of the SVI is employed.It is shown that, despite simpler construction compared to the 2-DVD, this instrument's data are suitable for such studies.
Even though this study is based on retrieval of ensemble mean snow density and not mass-dimensional relations directly, which could be more easily applied to radar retrievals and numerical weather prediction (NWP), there are a number of applications of such relations.Aikins et al. (2016) used ρ(D 0 ) to convert PSD observations to precipitation rate.Tong and Xue (2008), Dolan and Rutledge (2009), Matrosov et al. (2009), Huang et al. (2010) and Zhang et al. (2011) used mean snow density-median volume diameter relations for characterizing winter precipitation microphysics by radar.Kneifel et al. (2015) showed a connection between mean snow density and multi-frequency radar observations.Thompson et al. (2008) used the density relation by B07, and Iguchi et al. (2012) applied a similar density retrieval method to improve parametrization of snow microphysics in NWP models, for example.

Measurement setup
Measurements were made at the University of Helsinki Hyytiälä Forestry Field Station, Finland (61 • 50 37 N, 24 • 17 16 E), during the Biogenic Aerosols Effects on clouds and Climate (BAECC) field campaign (Petäjä et al., 2016) and during the consecutive winter of 2014/15.BAECC was a joint experiment between the University of Helsinki, the Finnish Meteorological Institute (FMI) and the United States Department of Energy Atmospheric Radiation Measurement (ARM) program.From 1 February through 12 September 2014 the second ARM Mobile Facility (AMF2) was deployed to the measurement site.The measurement setup was designed for snowfall intensive observation period of BAECC, called BAECC Snowfall Experiment (SNEX), which was undertaken from 1 February though 30 April 2014 and focused on measurements of snow microphysics.However, in order to extend the dataset, the measurements were continued upon completion of BAECC.In total, 23 snowfall cases from winters 2013/14 and 2014/15 were used in this study as summarized in Table 1.The snowfall cases were selected based on measurements of liquid water equivalent (LWE) precipitation accumulation by a weighing gauge, snow depth using a laser sensor and temperature measured by the automatic weather station of the FMI located 500 m distance from the measurement site.Only precipitation cases, where temperature was below or equal to 0 • C, were chosen, and the data were omitted when occasionally the temperature during the event rose above 0 • C. The experiments in both winters were organized in collaboration with the National Aeronautics and Space Administration (NASA) Global Precipitation Measurement (GPM) mission ground validation program.The surface precipitation measurements are carried out using a number of collocated instruments, such as NASA Particle Imaging Package (PIP), two OTT Pluvio 2 weighing gauges, two Parsivel 2 laser disdrometers (Tokay et al., 2014), a 2-DVD and a laser snow depth sensor by Jenoptik.To minimize effects of wind, a Double-Fence Intercomparison Reference (DFIR) wind protection (Rasmussen et al., 2012) was build on site as shown in Fig. 1 and discussed in more detail in Petäjä et al. (2016).Inside of the DFIR, the 2-DVD, one of the OTT Pluvio 2 s and one of the Parsivel 2 disdrometers were placed.In addition to the precipitation sensors, 3-D anemometers were deployed.The wind measurements were carried out at the heights of precipitation instrument sampling volumes.In this study data from the NASA PIP disdrometer and both OTT Pluvio 2 gauges are used.

Particle Imaging Package
The NASA PIP is the new generation of the SVI.The PIP, like the SVI, consists of a halogen lamp and a charge-coupled device full frame camera with sensor resolution of 640 × 480 pixels.The main differences between PIP and SVI are the camera and improved software.The camera is now capable of imaging with a frame rate of 380 frames per second, enabling measurements of particle fall velocities.The distance between the lamp and the camera lens is approximately 2 m.The lens focus is set at 1.3 m, where the field of view (FOV) is 64×48 mm, and the image resolution thereby 0.1×0.1 mm.The main advantage of PIP, as well of SVI, over other disdrometers is the open particle catch volume, which minimizes effect of wind on quantitative precipitation measurements (Newman et al., 2009).
The instrument records shadows of particles as they fall through the observation volume.Given the camera frame rate, multiple images of a particle are recorded and used to estimate its fall velocity.The depth of field (DOF) is determined by the processing software either rejecting or not detecting particles that are out of focus.Thus, the observation volume is defined by the FOV and the DOF.The expected particle size error due to the blurring effect is 18 % (Newman et al., 2009).From the recorded particle images a number of parameters describing particle geometrical properties are calculated with National Instruments IMAQ software.The measured diameter is given as the equivalent area diameter, which is the diameter of a circle with the same area as the area of a particle shadow.Other parameters, such as particle orientation, and bounding box width and height are also recorded.The aspect ratio of a particle is derived by fitting an ellipse to the bounding box utilizing the orientation of the particle.The aspect ratio is the minor axis in respect to major maxis of the fitted ellipse.The major axis also defines the minimum circumscribing disk, and the area ratio is defined as total area of shadowed pixels in respect to area of the circumscribing disk.

Weighing gauges and anemometers
The measurement setup includes two OTT Pluvio 2 weighing gauges, one inside and one outside the DFIR, with orifices of 200 and 400 cm 2 , respectively.There are differences in wind shielding as well.The Pluvio 2 200 is equipped with a Tretyakov wind shield and the Pluvio 2 400 with a combination of Tretyakov and Alter wind shields, as seen in the forefront in Fig. 1.
The gauges output several products of precipitation rate and accumulation.In this study, a non-real-time accumulation product is used as it is filtered for various sources of errors such as changes in the bucket mass due to evaporation, and as such should yield the most precise precipitation rate estimate among the output products.Because of the filtering, there is a 5 min delay in the recorded time series, which needs to be taken into account when comparing to other instruments.The precipitation accumulation values are recorded with a resolution of 0.001 mm, but non-real-time accumulation is output with a resolution of 0.05 mm.
The 3-D anemometer manufactured by Gill is located approximately at the height of the PIP on the field, respectively.The wind parameters, horizontal and vertical speed and horizontal direction, of Gill anemometer are measured every 10 s and averaged over 60 s.The mean and maximum of the 60 s wind speed averages and the mean wind direction for each event are given in the Table 1.

Snow depth sensor
The laser snow depth sensor, Jenoptik SHM30, is located on the measurement field, next to Pluvio 2 400.It is an optical sensor, which measures the snow depth by comparing signal phase information of the modulated visible laser light.It is a point measurement, and hence the piling of wind driven snow or random branches and leaves drifting on the snow pack can cause misreadings.To reduce this we have sheltered the measurement spot with a small wind fence and the instrument structure excluding the measurement pole is buried under the ground to prevent the piling of snow.The data are recorded every minute.

Retrievals of ensemble mean density, velocity-dimensional relations and PSD
Observations from the PIP and one of the weighing gauges are combined to retrieve snow ensemble mean density.Typically the gauge located inside of the DFIR, the Pluvio 2 200, is used for this retrieval.On a couple of days this gauge was not operational and data from the Pluvio 2 400 located outside of the DFIR were used instead.These dates are marked in Table 1 with asterisks in the LWE precipitation rate col- umn.As seen in the Table 1 the differences in accumulated LWE recorded by the two Pluvio 2 s are small, the largest being 15 %.Pluvio 2 200 inside the DFIR is typically measuring higher accumulations, which is expected because of the better wind protection.However, the observations do not show a clear indication that the observed precipitation accumulation difference depends on the wind speed.However, the difference seems to increase in respect to certain wind directions.
There are two openings from the measurement field, one to a road crossing (approx.130 • ) and the other to small field (approx.180 • ).If the wind is blowing from these directions the difference between the two gauges seem to increase.The retrieval procedure is described below and is similar to the one presented by B07, but with notable modifications.Prior to retrieval of ρ, PSD and velocity-dimensional relations are estimated.It was found, however, that the density retrieval is highly sensitive to the integration time.To minimize this, a variable integration time determined by the precipitation accumulation is used.The same integration time was applied to compute PSD parameters and v-D relations.

Particle size distribution
The PSDs are calculated from the PIP records of particles that fell through the observation volume.The observed distributions are defined with respect to equivalent area diameter D PIP , which is different from the apparent diameter of the 2-DVD and maximum particle dimensions used in other studies (e.g., Heymsfield et al., 2004).Wood et al. (2013) studied differences between diameter definitions and found that the diameter recorded by SVI is approximately 0.82 of maximum particle dimension.We performed a similar study by examining mean dimensions of rotated ellipsoids on a single projection, as shown in Fig. 2. The ellipsoids were defined by a long dimension a and a short dimension b lying nominally in the horizontal plane along the x and y-axes, respectively, and a short vertical dimension c lying nominally along the zaxis.The particle orientation was defined by Gaussian distri-bution of canting angles with a standard deviation of 9 • (Matrosov et al., 2005a) and a uniform distribution of azimuth angles.The equivalent area diameters D PIP of simulated particles were estimated from their projected areas onto the x-z plane and the resulting values were averaged over all orientations.The ratios of mean D PIP to the particle volume equivalent diameter, i.e., the diameter for which the particle volume V (D) = π 6 D 3 , for a number of combinations of vertical and horizontal aspect ratios are shown in Fig. 2. Assuming spheroids (Matrosov, 2007) and taking the typical vertical aspect ratio c/a = 0.6 (Korolev and Isaac, 2003;Matrosov et al., 2005b), we found that D PIP is roughly equal to 0.92 of a volume equivalent diameter.As can be seen, the conversion factor varies between 0.8 and 1.For ice particles with axis ratios smaller than 0.4, i.e., pristine ice crystals, this factor could approach 1.2.From this analysis we can conclude that the largest expected error is associated with observations of ice crystals.Dimensions of snowflake aggregates and graupel like particles are expected to be captured with a smaller error.In this study the same conversion factor of 0.92 is used for all the cases.As can be seen in Fig. 3 the median area and aspect ratios of the particles are 0.65 and 0.72, respectively.These observations also support our choice of a mean particle shape and the corresponding diameter transformation.Therefore, the results presented in the rest of the paper are using this volume equivalent diameter proxy.
Prior to calculations of PSD parameters, recorded PSD data are filtered to remove spurious observations of large particles.Following the procedure described in Leinonen et al. (2012), records of large particles were ignored if there was a gap of more than three consecutive PSD diameter bins.The bin size was set to 0.25 mm during the BAECC experiment and it was reduced to 0.2 mm for the winter 2014/2015.The PIP resolution is 0.1 mm and the minimum detectable particle diameter is approximately 0.3 mm (Newman et al., 2009).The smallest diameter bin used in calculations is 0.25 to 0.5 mm during BAECC and 0.2 to 0.4 mm in the following winter.
The PSD parameters were calculated using method of moments and assuming that PSD follows gamma functional form; see for example Ulbrich and Atlas (1998) and citations therein.The normalized gamma distribution N (D) in mm −1 m −3 was adopted following Testud et al. (2001), Bringi and Chandrasekar (2001) and Illingworth and Blackman (2002): with N w in mm −1 m −3 being the intercept parameter, D 0 the median volume diameter in mm, the slope parameter in www.atmos-meas-tech.net/9/4825/2016/Atmos.Meas.Tech., 9, 4825-4841, 2016 mm −1 and µ the shape parameter.Using the second, fourth and sixth moments for the non-truncated gamma PSD, M 2 , M 4 and M 6 , the PSD parameters were estimated as follows:

Ensemble mean density retrieval
The integration time, τ (t), of the ensemble mean density retrieval is driven by precipitation measurements of the Pluvio 2 .The step of the non-real-time accumulation output is 0.05 mm, causing the output interval to be on the order of several minutes even at moderate snow rates.With a short fixed integration time in timescales of minutes or tens of minutes, the produced ensemble mean density estimation would hence be more unstable, the lower the precipitation rate.Therefore, variable length time intervals driven by the gauge output are used with a selected threshold value of 0.1 mm.This corresponds to a τ (t) of 6 min for a LWE precipitation intensity of 1 mm h −1 .Effectively, the temporal resolution of the ensemble mean density retrieval is increased with increasing precipitation intensity, and in the analysis of the snowfall events in Table 1 the median τ (t) was 5 min.
As the integration time τ (t) is effectively driven by precipitation intensity, there is less variation in number of particles between intervals compared to a fixed time interval approach.With the selected accumulation threshold there are typically between 10 3 and 10 4 particles within a given integration time interval.However, with low precipitation intensities, τ (t) increases up to 1 h and retrieved ρ becomes less representative for the time interval in question.With LWE precipitation rates lower than 0.2 mm h −1 , the resolution of Pluvio 2 LWE measurements is insufficient and calculations of ρ become overly sensitive to recorded number concentrations.Correspondingly, similar unwanted sensitivity to LWE precipitation accumulation occurs when the number of particles observed by PIP within τ (t) is less than 800.Therefore, time intervals with precipitation rates or particle counts lower than these thresholds are excluded from our analysis.
Given a population of solid precipitation particles with volume equivalent diameters D over the integration time τ (t), the liquid equivalent precipitation accumulation in millimeters is approximately where ρ is the volume flux weighted population mean snow density in g cm −3 , ρ w = 1 g cm −3 is the density of liquid water, N (D, t) is mean particle number concentration over the integration time in mm −1 m −3 , v(D, t) is particle velocity relation in m s −1 and [D min , D max ] is the size range of snowflake observations from a disdrometer.From Eq. ( 8) we can estimate volume flux weighted snow density for each observation time interval as , using liquid equivalent precipitation accumulation G(t) as measured by the Pluvio 2 gauge, and retrieving averaged N (D, t) and volume flux using fitted v(D, t) based on measurements by the PIP as described in the following sections.
It should be noted that, unlike in the retrieval of PSD parameters where gamma PSD was assumed, ρ was retrieved without making any assumptions on the shape of the PSD distribution and instead measured PSDs are used in the calculations.

Comparison of derived mean density to snow depth observations
The definition of ensemble mean density here is the same as for mean bulk density in B07.They determine the densities for 5 min precipitation volumes derived with a 2-DVD disdrometer observations together with precipitation mass measured by a weighing gauge.B07 defined the volume of a single particle by summing coin-shaped sub-volumes together, estimated separately for both orthogonal projections and taking geometrical mean.As the diameter used in our study is the estimated volume-equivalent diameter, our results are comparable to B07.In Heymsfield et al. (2004), the volume of a single particle is defined as a function of circumscribing maximum diameter, and the population mean effective density is determined from ice water content.The estimated ensemble mean snow density is volume-weighted and expected to have lower values than the velocity-weighted snow density.The difference is not generally prominent especially with low-density aggregates, whose velocity-dimensional dependence is weak.
It should be noted that the derived density is inversely proportional to the snow ratio, R s , assuming that issues related to packing of snowflakes on the ground can be ignored.The snow ratio (Power et al., 1964;Ware et al., 2006) is used by operational weather services to estimate change in snow depth from LWE observations and can be defined as follows: where ρ(t) is the volume flux weighted snow density derived as shown in Eq. ( 9), P is the packing efficiency of snowflakes and C is the snow compression.Assuming that the packing and compression terms, or their product, are close to unity, the derived density can be tested against the commonly used assumption that 1 mm of LWE accumulation corresponds to 1 cm change in snow depth.In Fig. 4 the combined distribution of estimated snow ratios on temporal scales defined by the gauge accumulation for all the 23 events analyzed in this study is shown.It can be seen that the mean and median values, equal to 10 and 9, respectively, are very close to the commonly assumed value.This analysis assumes that packing efficiency of snowflakes is 100 % and compression of snow on the ground can be ignored or that snow compression counteracts reduction in snow density due to packing.The packing efficiency of snowflakes on the ground is not known.Random packing of the same size spheres has density of 64 % and dense packing of such spheres uses 74 % of the volume, corresponding to P = 0.64 and 0.74, respectively.Packing efficiency of equal spheroids depends on axis ratios, exceeding that of spheres, and could exceed 77 % (Donev et al., 2004).It is not unreasonable to expect that irregular shaped particles of variable sizes, such as snowflakes, would pack more efficiently than equal spheroids.At least, packing efficiency in excess of 90 % can be expected for spheres of several radii (de Laat et al., 2014).The packing efficiency of 70 % would mean that density of freshly fallen snow would be 30 % lower than that of falling snowflakes.The packing efficiency of 80 % would correspond to 20 % bias in estimated snowflake density from snow depth measurements or in 25 % underestimation of the snow depth change by using ρ(t).We do not know the exact value of the snow packing, but we could expect that in the worst case scenario it is about 70 % and probably closer to 80 % or even higher.It should also be noted that the snow compression would counteract this, but we are considering only freshly fallen snow and expect that the compression factor C is very close to unity.
One of the major uncertainties in the density retrieval is the assumption about particle volume.In this study we have assumed that snowflakes are spheroids with axis ratios of 0.6.Given this assumption, a conversion factor relating volume equivalent and observed disc equivalent diameters was defined.Figure 2 shows that for a reasonable range of ellipsoid axis ratios this conversion factor can range between 0.8 and 1.This range of values implies that the uncertainty in the density estimation can range from an overestimation by as much as 50 % to an underestimation by about 20 %.This range of uncertainty is much larger than what is expected from a comparison of the retrieved volume-flux weighted density and snow depth measurements, as was discussed previously.Therefore, by comparing the PIP derived and the directly measured snow depths, the validity of the derived values of ρ, and assumption of particle shape, can be checked.In Fig. 5 hourly change in the snow depth measured by the Jenoptik SHD30 is compared to the PIP derived snow depth.It can be seen that the agreement is good, with RMSE of 0.30 cm, linear correlation coefficient of 0.88 and normalized bias as low as −0.06.This comparison also gives confidence about the validity of the derived ensemble mean densities.

Effect of PSD truncation on derived ensemble mean snow density
The observed PSDs are truncated on left and right sides (Ulbrich and Atlas, 1998).They are truncated on the right side because of the instrument finite sampling volume and because natural sizes of hydrometeors do not extend to infinity.The truncation on the left, on the small-diameter side, is due to instrumental limitations and possible wind effects (Moisseev and Chandrasekar, 2007).Ulbrich and Atlas (1998) have presented a comprehensive analysis on how the right-side truncation affects the derived gamma PSD parameters.A similar study on the effects of the left-side truncation and other instrumental effects was presented by Moisseev and Chandrasekar (2007).Here we apply the method presented by Moisseev and Chandrasekar (2007) to estimate impact of PSD truncation on the derived mean snow density.
To investigate the impact of the PSD truncation on retrieval of mean snow density, a simulation study was performed.To initiate the simulation, the PSD parameters N w , D 0 and µ, together with parameters of m-D and v-D, are used.During the study it was found that the density estimation error is most sensitive to D 0 and µ and virtually independent of the other input parameters.Therefore, the results presented here assume that N w is constant and equal to 10 4 mm −1 m −3 .Further, only one m-D relation representative of all BAECC cases, as presented in Sect.4.3.1, is selected, and v-D representative of the snowfall with mean density ranging between 100 and 200 g cm −3 is utilized.The D 0 values were varied between 0.5 and 4 and µ values between −0.9 and 3.
At the first stage of the simulation, the number of observed particles was computed, assuming a Poisson distribution, with the expected number of particles being determined by PIP sampling volume and τ (t).Given this number of particles, their diameters were found by sampling a gamma probability density function, parameters of which are determined by the input PSD.To simulate the left-side truncation all particles with diameters smaller or equal to 0.25 mm, the PIP sensitivity threshold, were rejected.The rightside truncation was achieved by rejecting particles with sizes exceeding 3D 0 .For each D 0 and µ pair, 50 simulated PSDs were computed.Given the simulated truncated PSD the density is estimated in the same way as was presented above.This estimated density is compared to the one that is directly derived from the simulation input parameters and the results of their comparison is shown in Fig. 6.As one can see, the derived ensemble mean snow density is biased.The bias is largest for small D 0 , which is explained by the left-side PSD truncation.For D 0 larger than 1 mm, the bias decreases and approaches 2 %.Given that the error associated with PSD truncation is rather small for D 0 > 1 mm, and that most of the observations fall within this range, the truncation error is not corrected in this study.

Velocity-dimensional analysis
For the retrieval of volume flux weighted snow density, velocity-dimensional relations of falling snow need to be estimated.For each integration time interval, v(D) = a v D b v is fitted for velocity-diameter data from the PIP.The v(D) power-law fits to unfiltered data tend to be strongly biased by outliers.To address this problem, Gaussian kernel density estimation (KDE; Silverman, 1986) is used to find the most probable velocity for each diameter bin, and only observations with velocities within half width at half maximum from the bin peak KDE value are included in calculating the fit.Using the linear least squares method, a fit is performed for the data points in log-log scale to derive a power-law relation.It should be noted that using linear regression in log-log space does not optimally minimize residuals in linear space, but the method is used here as it does not overly emphasize the large end of the size spectrum.The retrieved velocity fits are shown for selected integration time intervals of the 18 March 2014 and the 22-23 January 2015 cases in the bottom of Figs. 7 and 8, respectively.
It should be noted that the power-law model, albeit widely used, may not necessary represent correctly velocities of ice particles over the complete range of diameters (Mitchell and Heymsfield, 2005).In many cases the fit can also be uncertain either because of narrow PSD or in presence of multiple particle types.During the 18 March, Finland was covered in a continental polar air mass.In the morning, a warm occluded front associated with a weak low pressure center approached Southern Finland from the southwest, bringing light snowfall.In the afternoon, Hyytiälä was in the warm sector of the frontal system and the relative humidity dropped, halting the snowfall around 12:00 UTC.Later in the evening there was a 1 h snowfall from a squall line, associated with a cold front passing over Southern Finland.
Time series of LWE snow rate, ensemble mean density and PSD parameters for the 18 March case are shown in Fig. 7.The bottom panels show measured fall velocities for selected integration time intervals, representing observations with different ensemble mean densities.Between the red dotted lines is the region where KDE is higher than half maximum for a given particle size.The fits are applied for data points between these lines.There is considerable scatter in particle fall velocity throughout the case and a bimodal PSD is present momentarily in the morning as can be seen in fall velocity panel Fig. 7a.
During the snow shower in the evening, liquid equivalent precipitation rates were recorded on average roughly 3 times more intense than earlier during the day, allowing retrievals of ρ and PSD parameters at high time resolutions.Strong short timescale variations of ρ and PSD parameters are recorded during this shower.The lowest ensemble mean density value of the case, 0.035 g cm −3 , is retrieved for time interval from 16:35 to 16:39, with concurrent D 0 value of 5.5 mm and N w of roughly 700 mm −1 m −3 .The corresponding fall velocity distribution visualized in Fig. 7b

22-23 January 2015
During 22 January 2015, similarly to the 18 March 2014 event, a warm occluded front associated with a weak low moved northwards over the Gulf of Finland.However, due to a blocking high over northwestern Russia, the low and the associated front were sustained over Southern Finland for the whole day of 23 January, causing weak continuous precipitation in the area.
Time series of LWE snow rate, ρ and PSD parameters for the 22-23 January 2015 case, with velocity-diameter fits from selected time intervals, are shown in Fig. 8.The case is characterized by continuous snowfall at LWE precipitation rates lower than 1 mm h −1 throughout the case.The velocity distribution for a given time interval has substantially less scatter compared to the 18 March 2014 case.The evolution 7. Evolution of snowfall intensity, ensemble mean density and particle size distribution parameters during 18 March 2015 with associated (v, D) from three selected time intervals (highlighted in gray).The red dashed lines mark the upper and lower velocity limits where for a given D the KDE value is higher than half maximum. of ρ and N w , as shown in Fig. 8, shows considerable similarities, suggesting a strong correlation.
The velocity-diameter fits shown represent a low ensemble mean density (ρ = 0.05 g cm −3 ) time interval of 01:03-01:16 (Fig. 8b) and two intervals of 22:30-22:52 and 02:06-02:14 (Fig. 8a, c) with higher values of ρ, 0.10 and 0.12 g cm −3 , respectively.Notable is the higher modal fall velocities and the absence of particles larger than 3 mm in the high density time intervals compared to the distribution in Fig. 8b.

v-D and density
In Fig. 9, particle fall velocity versus diameter data points combined from all the cases in Table 1 are divided into three categories according to the snow ensemble mean density of the time interval during which particles were observed.A least squares fit is applied to observations in each ρ range using the same procedure as for velocity-dimensional fits for integration time intervals, as described in Sect.3.5.The total number of observed particles is roughly 4 440 000, and for Figure 8. Evolution of snowfall intensity, ensemble mean density and particle size distribution parameters during the night between the 22 and 23 January 2015 with associated (v, D) from three selected time intervals (highlighted in grey).The red dashed lines mark the upper and lower velocity limits where for a given D, the KDE value is higher than half maximum.each density category numbers of particles included in the fitting process (within the red lines in Fig. 9) are approximately 1 140 000, 1 190 000 and 360 000, respectively.The fitted relations for ensemble mean density ranges are v(D) = 0.895D 0.244 0.1 < ρ ≤ 0.2 g cm −3 and ( 12) with RMSE values of 0.30, 0.30 and 0.35 m s −1 , respectively.
The coefficient is increased with density indicating higher fall velocities with more dense particles.There is also a clear increase in the slope of the fitted curve from the lowest density range to the 0.1-0.2g cm −3 range indicated by the increase in the power term.With particles in the highest density range the observed size distribution is narrow and hence the correlation between particle size and fall velocity is weak, and it is difficult to find an unambiguous relation between them.All things considered, the results are in line with the conclusion made by Barthazy and Schefold (2006) that the prefactor and power terms increase with riming degree, which in turn are strongly connected with density (Power et al., 1964).
Considering the definition of the volume equivalent diameter, relations in the form of Eqs. ( 11)-( 13) should be ideal for velocity-dimensional parametrization of radar observations as the average size of hydrometeors as observed by radar are largely defined by their volumes rather than their shapes.

Connection between PSD parameters and density
From the analysis of PSD parameters and their relations to ensemble mean density we have excluded data points representing integration time intervals where D 0 < 0.6 mm, as lower values of median volume diameter would imply that a substantial fraction of particles are too small to be observed with PIP.Applying this restriction, along with minimum thresholds set for particle count and LWE precipitation rate in density retrievals, as described in Sect.3.2, all in all 101 time intervals were discarded from the total of 1141 intervals of observations, 7173 min of snow observations for analysis.

Density and D 0
In Fig. 10, observed distributions of D 0 for the three different density regimes are shown.For the low-density particles, the maximum D 0 value rarely exceeds 5 mm, which is in agreement with observations of snow aggregates presented by Lo and Passarelli Jr. (1982).It can also be seen that D 0 distribution depends on density.Low-density particles are generally larger and vice versa.This dependence of D 0 on ensemble mean density is not surprising, given that they are related as was previously shown by B07 and discussed in more detail below.
Relation between ρ and size (D 0 ) is illustrated in Fig. 11.The areas of individual data points are proportional to the particle counts of the corresponding observation time intervals.The overlaid black solid curve, a least squares fit applied for all cases in Table 1, is given by where D 0 is in mm and ρ is in g cm −3 .As the two examined winters were seen to have notable differences between each other in the snowfall type and average ρ, corresponding relations were also calculated separately for the winters and are given by  easy to compare.Especially Eq. ( 16) is in good agreement with B07's results.The ensemble mean density is on average higher for snow events recorded during BAECC, which suggests more riming occurred during those events.Indication to this is that the ARM AMF2 dual-channel microwave radiometer located on the same measurement field detected the presence of liquid water more than 80 % of the BAECC SNEX campaign time (Petäjä et al., 2016) and the presence of supercooled liquid layers could also be observed in the backscatter coefficient and circular depolarization ratio measurements of the co-located ARM HSRL (High Spectral Resolution Lidar) in the majority of the BAECC cases (Goldsmith et al., 2014).In general the BAECC winter was milder than the next winter 2014-2015, and the case duration weighted average of maximum recorded temperatures was almost 1 • C higher for BAECC events compared to the value for winter 2014-2015 cases.The temperatures closer to 0 • C could mean increased aggregation as stated in B07, and therefore decreased density values, as well as different snow habits compared to colder cases.The mass-dimensional relation in power-law format m = a m D b m can be derived from the retrieved ρ-D 0 relations (Eqs.14-16) by assuming gamma PSD and describing the ensemble mean density approximately as The integration limits are defined from 0 to infinity for deriving the analytic solution, though the true range is narrower because of left and right truncation of the observed size spectrum.As shown in Fig. 6, the ensemble mean density is overestimated because of the truncation.The estimation bias of density ranges between 20 % for D 0 smaller than 0.75 mm and about 2 % for D 0 larger than 2 mm.Since for the estimation of the m-D relation most of the observed D 0 values are higher than approx. 1 mm as shown in Fig. 10, there is only minor contribution of the smaller D 0 values, and we assume our error in ensemble mean density to be close to 2 % W Figure 11.(D 0 , ρ) for all cases listed in Table 1.Area of each dot is proportional to the number of particles in corresponding integration time interval.Power-law fits are shown separately for BAECC winter cases (blue) and cases from the following winter (green).because of truncation.This corresponds to an error of 2 % in the prefactor a m as well, if it is assumed that the truncation does not introduce significant changes in the exponents of the ρ-D 0 and m-D relations.
Taking the three velocity exponents from Eqs. ( 11)-( 13), and assuming exponential PSD, the derived prefactors and exponents of mass relation are shown in Table 2, having the volume-equivalent diameter proxy in millimeters and mass given in grams.The factor 0.1 in Eq. ( 18) is derived from unit conversion, as ρ is in g cm −3 .The values of prefactor a m are not sensitive to the changes in the velocity exponent b v (changes in b v are resulting less than 1 % deviation a m values), though there is a small increase in a m with increasing b v .The prefactor is more sensitive to shape parameter µ of the gamma PSD; the value of a m increases by 24 % as µ is increased from 0 to 1.With value of µ = 3 the increase in the prefactor a m value is 48 %.The shape factor of snow PSD is 2.04 7.5814 Mitchell et al. (1990) 2.0 3.0926 Locatelli and Hobbs (1974) 1.9 5.1134 known to be noisy and thus often exponential distribution is assumed.With b v = 0.217 the derived mass-dimensional relations for all cases and for both studied winters separately are plotted against literature values in Fig. 12.The derived exponent b m for the studied cases is in line with literature values, close to 2, but the prefactor a m values are higher than in the presented relations in Table 3.The highest value of a m is for the BAECC cases, indicating conditions of riming.
The high prefactor values might manifest the Finnish winter conditions.Because of the vicinity of the Baltic Sea, the air is more moist than, for example, in continental conditions.

N w and density
Distributions of observed N w values also exhibit dependence of N w on the ensemble mean density, as shown in Fig. 13; i.e., N w increases with density.The modal values of N w are approximately 5000, 40 000 and 80 000 mm −1 m −3 for ensemble mean density ranges 0.0-0.1,0.1-0.2 and > 0.2 g cm −3 , respectively, with the vast majority of N w values spanning less than 2 orders of magnitude for a given ρ range.This dependence of N w on density is somewhat unex- pected.There is no obvious reason to expect that N w would depend on density.However, because D 0 and density are related, dependence of N w on density potentially arises from the dependence of N w on D 0 .
To verify this, the partial correlation analysis of the relation between log values of N w and density while controlling for log value of D 0 was carried out.It was found that there is a moderate negative partial correlation, −0.33, between N w and density while controlling for D 0 .However, the zeroorder correlation between N w and density is 0.52.The analysis confirms that the observed relation between N w and density is due to their relation to D 0 .It is not clear, however, what the meaning is of the found negative partial correlation between N w and density.
A relation between N w and snow particle size is shown in Fig. 14a.A linear least squares fit is applied for (D 0 , log(N w )), and the corresponding relation between N w and D 0 is given by N w = 2.492 × 10 5 × 10 −0.620D 0 .
(20) Bringi and Chandrasekar (2001) show that there is a weak tendency for N w to decrease with increasing D 0 for rain (their Fig. 7.17), but to our knowledge this is the first attempt to find a climatological relation between D 0 and N w for snow.It should be noted, however, that the observed relation is partially caused by data filtering, which removes low precipitation rate data.There is a high amount of scatter when N w < 1 × 10 3 mm −1 m −3 .The data points in this area are more contained when D 0 is multiplied with ρ   However, the difference in correlation coefficients for the fits in Fig. 14a and b, given by −0.87 and −0.85, respectively, is minimal.The lower scatter in Fig. 14b for N w in the sub 10 3 mm −1 m −3 range seems to be compensated by slightly more scatter in the higher end of the distribution.

PSD shape parameter, µ
In Fig. 15 the normalized frequencies of the gamma PSD shape factor µ are visualized in the three ensemble mean density ranges.Unlike D 0 and N w , µ does not seem to have a clear correlation with ensemble mean snow density, although a weak tendency for µ to increase with density is possible.Instead, the values of µ are scattered around approximately 0, with deviation increasing with density.In the ensemble mean density ranges 0.0 to 0.1 and 0.1 to 0.2 g cm −3 the kernel densities peak at −0.15 and 0.62, with standard deviations of 0.97 and 1.58, respectively.For the integration intervals with ρ > 0.2 g cm −3 , the distribution of µ is more spread, with standard deviation of 2.0 and median of 0.76.The observations support the findings of B07 and Heymsfield et al. (2008), who have found that low-density particles generally have exponential or slightly super-exponential distributions.This suggests that the exponential PSD would be most appropriate for describing low-density aggregated snow and less so when strong riming occurs.

Conclusions
Microphysical properties of snow in Southern Finland were documented using observations from PIP and a weighing gauge.The data were collected during US DOE ARM funded BAECC campaign and the consecutive winter.It is shown that there is a detectable difference in measured snow properties between consecutive winters.Snow observed during BAECC is denser than during the next winter.The derived m-D relations from two winters are also different, and the difference is namely in the prefactor of the power-law relations.
It is found that D 0 and N w parameters of gamma PSD are correlated with ρ.While the relation between ensemble mean density and D 0 is not surprising, since these two parameters are related, the correlation between N w and ρ is interesting.This correlation arises from the observed connection between N w and D 0 .It should be noted that this observed connection is partially due to data filtering that removes low precipitation rate data from the analysis.However, it indicates that for heavier precipitation aggregation is an important snow growth process.During snow growth by aggregation, N w should decrease while D 0 increases, as was found by (Lo and Passarelli Jr., 1982).The shape parameter of the gamma PSD, µ, does not seem to depend on ensemble mean density and its average value is close to 0, which is in line with studies reported in literature.
Dependence of v-D relation on ensemble mean density was also studied.It was found that the prefactor of the v-D power law depends on density.It is higher for higher densities.This result is in agreement with the conclusion made by Barthazy and Schefold (2006): the coefficient and power terms increase with riming degree.
The presented study uses the newly developed instrument Particle Imaging Package, which is a new generation of SVI.It is shown that data collected by this instrument are adequate for such studies.While the instrument only observes particle shapes projected to single 2-D plane, as opposed to 2-DVD or MASC, it has a larger sampling volume and its observations are less affected by wind (Newman et al., 2009).Additionally, the instrument itself is operationally more robust and requires less maintenance, enabling deployment in sites with remote locations and harsh field conditions.

Data availability
The data of the video distrometer (PIP), the precipitation gauges and the snow depth sensor used in this study are available at http://www.arm.gov/campaigns/amf2014baecc#data or by request from D. Moisseev (dmitri.moisseev@helsinki.fi).

Figure 1 .
Figure 1.Snow precipitation instruments on the measurement field in Hyytiälä.The Pluvio 2 200 is inside the wind protection on a platform and the PIP lamp can be seen at right on the ground.The view of the picture is to southwest and the distance from the platform to the treeline behind is approximately 20 m.

Figure 2 .
Figure 2. The ratio of the diameter observed by PIP, D PIP , to volume equivalent diameter D.

Figure 3 .
Figure 3.The distributions of snowflake (a) aspect ratio and (b) area ratio as observed using PIP with interquartile ranges visualized and median values shown.

Figure 4 .
Figure 4. Distribution of snow ratios, ratio of snow depth change to LWE, calculated from retrieved ensemble mean densities with interquartile range, and median and mean values.

Figure 5 .
Figure 5. Scatterplot of the hourly change of snow depth measured with Jenoptik SMH30 and estimated from volume flux using PSD and fall velocities as measured by PIP.The data include all the studied cases except 10-11 January 2015.

Figure 6 .
Figure 6.Computed normalized bias and standard deviation of estimated mean snow density as a function of µ and D 0 .The shaded area indicates data that are not included in the analysis because derived D 0 is smaller than 0.6 mm.The increased values of bias at low D 0 values is due to left side truncation of the observed PSD, which is caused by the instrument sensitivity.At larger D 0 values the bias approaches value of 0.02.
is characterized by low values of velocity fit coefficients a v and b v .Within the following 20 min, D 0 decreases down to roughly 2 mm, N w increases to 2 × 10 4 mm −1 m −3 and retrieved values of ρ peak at over 0.2 g cm −3 between 16:54 and 16:58 and again from 17:05 to 17:08.Corresponding fall velocity distribution between 16:54 and 16:56, shown in Fig. 7c, is characterized by substantially higher values of a v and b v than 20 min earlier, which possibly indicates the onset of riming.

Figure 9 .
Figure 9. Probability densities of (D, v) in three ensemble mean density ranges ([ρ] = g cm −3 ).Dashed lines mark the full width at half maximum KDE in each diameter bin.Power-law functions are fitted for data between those lines.

Figure 12 .
Figure 12.Derived m-D relations assuming exponential PSD in comparison relations presented in literature are shown in Table3.A conversion of maximum dimension to volume equivalent diameter is done by assuming axis ratio of 0.6.

Table 1 .
Liquid water equivalent precipitation accumulation measured with Pluvio 2 200 and 400, change in snow depth and maximum and minimum temperature, maximum and minimum relative humidity, mean and maximum wind speed and mean wind direction of the studied snow events.Events before the horizontal line are recorded during the BAECC campaign.

Table 2 .
The prefactors and exponents of m = a m D b m derived for exponential PSD with different values of exponent b v of velocity relation.The mass given in grams and the volume-equivalent diameter proxy in millimeters.