Laser attenuation in falling snow correlated with measurements of snow particle size distribution

Abstract. We investigate transmission through falling snow using a 3D imaging lidar at 1.56-μm wavelength. The lidar is based on the time-correlated single-photon counting technique. Experimental transmission data are compared with Mie theory transmission calculations based on snow particle size distribution simultaneously measured with a laser disdrometer. The calculations were performed in two ways, using the Beer–Lambert approach where all radiation interacting with a hydrometeor is considered extinct and an approach that includes effects of the forward scattering. Comparison of these two methods shows that inclusion of the contribution from forward scattering gives better agreement between experiment and calculations than using extinction only. When comparing the results using scattering calculations with a curve fit approach based on precipitation rate, it is evident that both the Beer–Lambert approach and the forward scattering approach give a much better correlation between experiment and calculations than relying on precipitation rate, as measured with the disdrometer, only.


Introduction
Optical attenuation in degraded visual environments, such as precipitation, can severely affect the performance of laser-based systems. Examples of technologies where a thorough understanding of signal behavior in precipitation is needed to predict the performance are lidar and free space optical communication (FSO), and considerable interest has been devoted to investigating precipitation effects here. Grabner and Kvicera 1 presented a multiple scattering model based on Monte Carlo simulations for calculation of signal attenuation by rain and fog for an FSO link. Further investigating signal attenuation within the FSO context, Korai et al. 2 presented work on a model for attenuation effects due to rain. In their work, they integrated FSO attenuation calculations with a global rainfall model, which given local rainfall statistics can generate precipitation rate preserving maps of rainfall events based on a number of synthetic rain fields. Goodin et al. 3 reported on a model for prediction of lidar performance degradation due to rainfall using analytical expressions for the drop distributions, and integrated this model into a 3D simulator for autonomous vehicles. Roy et al. discussed the interaction between falling snow and a 3D lidar scanner. In their work, they calculated the lidar signal in model snowfalls, and based on their findings they discussed a filtering algorithm for suppression of noise due to scattering from the snow particles. 4 Transmission of optical radiation through falling snow has been measured and discussed in a number of publications. O'Brien 5 reported on attenuation of visible light using both visibility observations and photometer measurements. Snowfall characteristics were obtained by visual inspection of snow particles and mass accumulation measurements. Miller investigated transmission of radiation from both HeNe and a CO 2 lasers through an artificially generated snow field. 6 The snow particle characteristics were measured with a digital camera and a PSD *Address all correspondence to Mattias Rahm, mattias.rahm@foi.se generated by post processing the images under the assumption of particles being spherical. Using the spherical particle assumption, the transmission was calculated with a geometrical optics model that included contribution from diffraction. Miller's result showed a spectral dependence of transmission, which, based on the model result, was attributed to the wavelength dependent diffraction lobe. 6 Mason reviewed experimental attenuation data from literature. 7 Mason's work showed that there is a wide variation in attenuation as a function of precipitation rate between different snowfalls. In the analyzed data, the extinction coefficient varied with a factor of five for a fixed precipitation rate, and the extinction coefficient for snowfalls could be more than one order of magnitude greater than the corresponding extinction coefficient in rain. The variation between different snowfalls was attributed to both variation of the particle size distribution (PSD) due to difference in meteorological conditions and to the presence of fog. Seagraves and Ebersole built on Masons review and presented an empirical model trying to relate attenuation in falling snow to meteorological conditions using additional data from a measurement campaign in Vermont. 8 Hutt et al. 9 discussed the spectral dependence of extinction measurements in snowfalls, both in the single and multiple scattering regimes. Their work indicated that the same spectral dependence of transmission reported by Miller is also present at longer propagation distances and in situations with significant multiple scattering effects. Nebuloni and Capsoni reported on measurements of laser signal attenuation due to falling snow in a research FSO link. 10 Common for these investigations are the noted large variation in extinction coefficients with respect to precipitation rate. The reason for the large variation lies in the complex processes of ice crystal growth and subsequent aggregation into larger structures, which are governed by the ambient meteorological conditions. 11 These variations in growth will be manifested by variations in the mass density and PSD.
The development and refinement of equipment for precipitation measurements, e.g., disdrometers that perform in-situ measurements of hydrometeor size and velocity, [12][13][14][15] have opened up possibilities for detailed correlation of precipitation characteristics with laser transmission measurements. Examples of the use of knowledge of the PSD when investigating transmission in precipitation are Brazda and Fiser, 16 and Peckhaus et al. 17 Brazda and Fiser investigated the impact of rain on an 850-nm FSO-link, and Peckhaus et al. investigated the effect of rain and falling snow on a high power laser beam centered around 1 μm.
In this paper, we report on an investigation of transmission through falling snow using a 3D lidar imager with a laser wavelength of 1.56 μm. Using an imaging lidar setup single path transmission cannot be experimentally separated. This setup, however, conveniently allows investigation of target areas larger than what is typically used in single path setups. The size of the target area can affect the atmospheric transmission by contribution from forward scattered radiation. The 3D imager used in the experiments relies on the time-correlated single-photon counting technique (TCSPC), and the snowfalls were characterized using a laser disdrometer. The experimental transmission data are compared and correlated with numerical calculations based on the detailed PSD measured by the disdrometer. The paper is organized as follows: In Sec. 2, the experimental setup and methodology is described. In Sec. 3, the equations employed to estimate transmission based on disdrometer PSD, using both extinction only and extinction plus forward scattering are derived. In Sec. 4, the results of both experiments and calculations are presented and compared. Sections 5 and 6 present discussion of the results and a summary, respectively.

Experimental Setup and Signal Processing
The setup used for transmission measurements is based around the time-correlated single-photon counting (TCSPC) 3D imaging system 3D-FALKen v3.11, shown in Fig. 1(a). The version 3.0 of the system is discussed in more detail in Henriksson et al. 18 The main components are the same in v3.0 and v3.11. The difference is that in the latter, the system is internally triggered from the camera, and neutral density (ND) filters can be used in front of the zoom lens. The TCSPC system consists of a Princeton Lightwave Inc. Falcon photon counting camera and a pulsed 1557-nm laser source. The camera has a 128 × 32 array of single-photon-sensitive Geiger-mode avalanche photodiodes (GmAPD) with 50-μm pixel pitch and readout electronics to measure the time of detection relative the camera trigger pulse. The array is mounted on a Kowa Z20750AMP zoom lens equipped with a 5-nm bandpass filter centered at 1558 nm for background suppression, an external motorized iris for control of effective aperture size, and a variable number of ND filters for additional control of the overall detector signal levels. The used focal length of 750 mm together with the pixel pitch give rise to an angular resolution of 67 μrad. The laser source is a BKtel HFL-310am-COL50 at 1557 nm with a pulse length of 0.9 ns, a variable pulse energy up to 14 μJ, and a variable pulse repetition frequency between 30 and 1000 kHz. The divergence of the laser beam is matched to the FOVof the camera equipped with a 750-mm focal length lens. Even though the initially linearly polarized radiation from the laser source undergo polarization changes when scattered from the snow particles and the target, no polarization filtering is performed on the detected signal. Consequently, with respect to the atmosphere, only absorption and spatial redirection of radiation due to scattering will contribute to the measured variation in transmission.
In brief, the TCSPC 3D imaging system measures the time between the trigger signal and the first absorption of a photon in each pixel. The temporal resolution of the 3D-FALKen is 500 ps. If consecutive pulses are sent out, a histogram over time from trigger pulse to photon absorption can be constructed. The histogram gives depth information in the pixel FOV, where the height and position of a peak in the histogram reveal reflection strength and distance to the object, respectively. An example of a histogram captured during the measurements can be seen in Fig. 2. Using several GmAPD in an array, a 3D image of a scene can be constructed.
The target used for transmission measurements was a 1.2 × 1.2 m 2 resolution test chart [shown in Fig. 1(b)] positioned at a distance of 708 m from the 3D-FALKen. In each measurement, the 3D-FALKen was scanned horizontally across the target four times with a frame rate of 80 kHz and a rotation rate of 1 deg ∕s. The data were subsequently assembled into panorama pixels and corrected for saturation effects. Panorama image generation is discussed in more detail by Henriksson et al. 18 and Henriksson and Jonsson,19 and the saturation correction by Jonsson et al. 20 A 2D panorama image of the target is shown in Fig. 1(c). The 2D panorama image was constructed by filtering each pixel histogram for peak value. From the panorama image a region of interest (ROI) covering 24 × 24 panorama pixels [indicated by the red box in Fig. 1(c)], which at the target corresponds to 1.1 × 1.1 m 2 , was used for the transmission measurements. The ROI was chosen to cover as much of the target as possible without including effects from the target edges. A number of pixel rows in the 3D-FALKen exhibit regions with non-uniform pixel sensitivity. These rows were removed from the data in order not to affect the measured signal. The position of the target in the panorama image varied slightly between measurements. Consequently, the effective number of pixels in the ROI varied between 16 × 24 and 19 × 24 in the measurements. The frame rate, the rotation rate, and the number of scans resulted in a 3D panorama image of the ROI containing around 39,000 measurements per panorama pixel.
In each pixel histogram the background signal was calculated as the mean value from start to 10 ns before the target peak (an example is shown in Fig. 2), corresponding to a background integration time of around 260 ns. The background signal consists of detector dark counts, ambient light reflected by the target, and ambient light and laser radiation scattered by aerosols and snow particles inside the beam path. The background was subtracted from the data and a target signal was calculated as the sum of counts in a temporal region of AE2.5 ns around the peak, i.e., 11 time bins in the histogram. A mean signal was obtained by calculating the average target signal inside the ROI for all four scans. The reflection from the whole target was used when calculating the signal. Therefore, the black bars used for resolution estimation, evident in Fig. 1(b), will result in an overall reduction of the reflection cross-section, affecting all measurements in the same way.
To relate the measured signal to attenuation induced by precipitation in a geometry dictated by a lidar system, the starting point is the lidar equation. The lidar equation, which is described in detail elsewhere, 21 can in the present experimental geometry be written E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 3 0 7 In Eq. (1) N D is the number of detected photons, N E is the number of emitted photons, σ is the target reflection cross-section, A illum is the illuminated area, A rec is the effective aperture area of the detector, L is the distance from lidar system to target, η sys is the receiver system efficiency, and η atm is the one-way atmospheric transmission efficiency. The laser and imager were scanned across the target when creating the panorama image. Thus, depending on the laser pointing direction, the finite width of the laser beam will result in a varying illumination of the target. In all measurements, the imaged region stretched well beyond the target in the horizontal direction. Therefore, the varying illumination within a scan, which will enter Eq. (1) as a reduction in N E , will act as a constant modification to the equation. It can thus be neglected when comparing repeated measurements of the target. The varying illumination in the vertical direction can affect signal levels through changes in panorama image center with respect to the target. However, in our measurements, the signal variation due to vertical movement between the different measurements was insignificant compared to other signal variation contributors and is therefore not included in the analysis. The atmospheric transmission efficiency is affected by scattering and absorption in the atmosphere, which in the present case involve radiation interaction with aerosols and hydrometeors. For the used wavelength, the atmospheric gas absorption is negligible. When comparing repeated measurements of the target σ, A illum , and L will remain constant. Depending on the ambient conditions, different number of ND-filters were placed in front of the lens. These variations will affect η sys , but other aspects of η sys , e.g., detector quantum efficiency and lens losses, will remain constant. Also the iris diameter, affecting A rec , was changed depending on ambient conditions. Therefore, a scaled version of Eq. (1), the scaled signal S scaled , can be introduced as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 6 7 5 In Eq. (2), S scaled is the number of detected photons scaled with variations in laser pulse energy, variations in effective aperture area, and variation of system efficiency through changes in optical density (OD) in the ND filters. C is a constant including contributions from the target reflection cross-section, the illuminated area, the distance between lidar and target, and the part of the system efficiency being constant during the measurements. This scaled signal directly relates to the atmospheric efficiency and will be used in the comparisons. The precipitation was characterized with a Thies 5.4110.01.200 laser disdrometer positioned close to the 3D-FALKen. The disdrometer simultaneously measures hydrometeor diameter and fall speed from which, together with temperature, precipitation type can be deduced. Data are binned into 20 velocity classes (0 to >10 m∕s) and 22 diameter classes (0.125 to >8 mm) with a temporal resolution of 1 min. The operating principle is a laser plane being attenuated by the falling hydrometeors. By the strength and duration of the attenuation, the size and velocity of the hydrometeor can be inferred. Precipitation rate (mm/h) is also reported by the disdrometer. Details of the algorithms classifying the precipitation type, inferring diameter, and precipitation rate are unfortunately not reported by the manufacturer. To use the precipitation characterization in calculations of transmission, the size/velocity histogram can be converted to a PSD. If M v;i denotes the number of hydrometeors in velocity class v and diameter class i, the per-diameterclass volumetric number density, N i (m −3 mm −1 ), can be calculated as 22 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 4 1 4 In Eq. (3), A i is the effective sampling area (m 2 ) of diameter class i, ΔD i is the width of the i'th diameter class (mm), Δt is the integration time (s), and V v is the velocity of the v'th class (m/s). Recently, Fehlmann et al. evaluated a Thies laser disdrometer similar to the one used in the reported campaign with respect to monitoring rain and snowfall events. 15 In the experiments, they compared the results from the laser disdrometer with the results from a pluviometer from OTT and a 2DVD video disdrometer from Joanneum research. The pluviometer was used as a reference for precipitation rate and the 2DVD disdrometer was the reference for PSD measurements. They could conclude that in their snowfall experiments the laser disdrometer slightly underestimated the number density in the PSDs compared to the 2DVD. With respect to recorded precipitation rates, Fehlmann et al. could see no clear bias with respect to the pluviometer, but instead a higher degree of uncertainty in the precipitation rate estimation from the laser disdrometer compared to the reference equipment. The higher uncertainty was attributed to underlying assumptions in the laser disdrometer built-in precipitation rate estimation algorithms. It has been noted in literature that snowfall sometimes occur in the presence of fog. 5,7 Therefore, to characterize the atmosphere with respect to particles smaller than the resolution limit of the disdrometer, an Alphasense OPC-N3 optical particle counter was used. The OPC-N3 measure the scattered light by individual particles carried in a sample air stream through a laser beam. The OPC-N3 does not dry the particles prior to measurement, and this feature makes the device suitable for size measurements of fog droplets or aerosols with adsorbed water. The results are binned into 24 diameter classes ranging from 0.35 to 40 μm, and data are recorded on a per second basis. During the experiments, to get reliable sampling of the particle content, a 10-min rolling mean value was used when calculating the aerosol PSD. It can be noted that there remains an interval between 40-and 125-μm diameter that is not covered by the sensors. However, the presence of fog or lack thereof can be determined with this particle counter.
In addition to the laser disdrometer and the optical particle counter, the visibility was simultaneously measured with a forward scatter visibility sensor from a Vaisala PWD22. The sensor measures visibility in the range of 0 to 20 km.

Numerical Calculations of Precipitation Induced Signal Attenuation
The starting point for evaluation of transmission efficiency of laser propagation in the atmosphere is often the Beer-Lambert law, which can be understood in a geometry depicted by Fig. 3. The Beer-Lambert law, assuming a homogenous attenuating atmosphere, states that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 6 2 5 Here EðxÞ is the energy of the pulse at a position x, E 0 is the initial energy of the pulse, and B ext is the extinction coefficient, which is the sum of the absorption and scattering coefficients.
In this work, B ext is calculated using Mie theory described below. In the Beer-Lambert law, all radiation that have interacted with the atmosphere via absorption and scattering are considered extinct. The atmospheric transmission efficiency, i.e., the ratio of energy at the detector and at the source, can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 5 2 3 In Eq. (5), the subscript BL has been added to the transmission efficiency to emphasize that it is calculated in accordance with the Beer-Lambert law. However, as shown by, e.g., Miller 6 and Hutt et al. 9 including forward scattering within the solid angle subtended by the detector aperture in the expressions for transmission in precipitation will in many cases give better agreement between calculations and experiments than using only the Beer-Lambert law. If the scattering is assumed to take place within the single scattering regime, Eq. (5) can easily be modified to approximate effects of forward scattering. The single scattering regime implies that a photon can only be scattered once. Photons that have previously been scattered within the detector solid angle will not participate in any further scattering events. Consequently, by adding an expression of the integrated radiation scattered within the detector solid angle, the expression for atmospheric efficiency can be written E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 3 5 1 Ω d ðxÞ represents the solid angle subtended by the detector aperture at position x and β sca is the angle dependent differential scattering cross-section. Furthermore, in Eq. (6), the subscript FS has been added to the atmospheric efficiency to emphasize that it is including effects of forward scattering. It is evident that within the single scattering regime, including forward scattering in the atmospheric efficiency calculations will always increase the transmission. Snow particles are not spherical in shape. The used disdrometer, however, report hydrometeor sizes in terms of diameter. Therefore, as a first-order approximation, motivated by the random orientation of falling snow particles with respect to the incident radiation, the hydrometeors are treated as spheres characterized by their diameter. Treating the hydrometeors as spherical objects allows Mie theory to be used in calculations of scattering and extinction parameters. The Mie calculations were performed with the freely available code package developed by Schäfer, 23,24 and refractive index data for ice was taken from Warren. 25 Based on the size dependent extinction and scattering cross-sections calculated using Mie theory, and ignoring the negligible gas absorption at this laser wavelength, the extinction coefficient and the differential scattering cross-section for the measured PSDs can be calculated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 6 9 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 4 4 In Eq. (7), r represents the particle radius, σ ext ðrÞ the size dependent extinction cross-section for a single particle calculated with the Mie code, and nðrÞ the measured PSD. In Eq. (8), σ sca ðθ; rÞ is the size and angle dependent differential scattering cross-section for a single particle, calculated with the Mie code. Using the rotational symmetry of the differential scattering cross-section imposed by the assumption of spherical particle shape, the expression for the transmission efficiency including forward scattering contribution can be written E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 5 5 0 sinðθÞβ sca ðθÞdθ dx: In Eq. (9), θ d ðxÞ is the half angle subtended by the detector aperture at a position x, shown in Fig. 3. Here, a zero divergence laser beam aimed at the center of a circular aperture is assumed. This is of course an approximation of the measurements where the divergent laser beam is swept across the square target to produce the panorama image. For the case of scattering toward the circular detector aperture, the integration is straight forward. For the case of scattering toward the square target, the square is approximated by a disk with equal geometric cross-section. The effect of the difference in geometry is expected to be small as the difference in angles is low for most of the path. Only when x is very close to L the difference in angles becomes substantial. Furthermore, effects of the spherical nature of the diffusely reflected waves from the target toward the circular detector aperture, and the detector's limited field acceptance angle, have both been neglected. Including the contribution from forward scattering in the expression for transmission efficiency breaks the symmetry of Eq. (1). Transmission efficiency of the outbound radiation path will not necessarily be the same as the transmission efficiency of the inbound path. The reason is differences in target and aperture sizes. If the target is larger than the aperture of the detector, a larger fraction of the radiation scattered off-axis will make it to the end of the outbound path than to the end of the inbound path, resulting in a higher transmission in the former case. Thus, Eq. (2) must now be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 2 9 0 S scaled ¼ Cη atmFS−out η atmFS−in ; where η atmFS−out and η atmFS−in are the outbound and inbound transmission efficiencies, respectively. To compare calculations with experimental data, the constant C, representing effects of target reflection characteristics and the constant part of the system efficiency (e.g., lens losses and quantum efficiency), must be determined. In this work, C has been calculated as the median value of C ¼ S scaled ∕η 2 atm for the experimental data points in the case of Beer-Lambert law type transmission, and the corresponding expression based on Eq. (10) in the forward scattering case.

Results
The measurements were performed in Älvdalen, Sweden, during March 3, 2020, and March 4, 2020. The weather was characterized by temperatures in the range of −4°C to −1°C and snow showers with varying precipitation rate, as shown in Fig. 4. During the campaign, transmission was measured with the imaging lidar at 16 different times, spanning recorded precipitation rates from 0.00 to 2.44 mm∕h, also indicated in Fig. 4. Figure 5 shows a scatterplot with simultaneous measurements from the visibility sensor and the disdrometer. Gray markers show all recorded data from the campaign, and red markers indicate data points where transmission was measured. As can be seen, there is considerable variation in visibility with respect to the recorded precipitation rates. To investigate the influence of the aerosol concentration, transmission calculations were performed using data from both the disdrometer and the particle counter. The calculations show that for the cases when transmission was measured with the imaging lidar, the aerosol concentration (e.g., due to fog) was too low to significantly affect the transmission efficiency results. Therefore, the variation in measured transmission can be attributed to precipitation only. Figure 6 shows the calculated transmission efficiencies based on the disdrometer PSDs, using both Eqs. (5) and (9). As expected, the Beer-Lambert approach (η atmBL ) results in the lowest transmission since all radiation interacting with the snow is removed from the signal. The approach, including radiation scattered within the target/detector solid angle (η atmFS ), shows significant difference between outbound and inbound transmission. The difference between outbound and inbound transmission is due to difference in size between the target and the 3D-FALKen aperture, where the ROI of the target has an area of around 1.2 m 2 and the 3D-FALKen iris has an area of around 7 × 10 −3 m 2 when fully open. The solid angle subtended  by the target allows a larger fraction of the radiation scattered in the forward direction to contribute to the signal than on the way back toward the detector. Comparing the Beer-Lambert result to the result from the approach including forward scattering, one can see that the inbound transmission in the forward scattering case closely resembles the Beer-Lambert result. The solid angle subtended by the 3D-FALKen iris is small enough to block the majority of the forwardly scattered radiation. The outbound transmission also closely resembles the Beer-Lambert result at very low precipitation rates. At around 0.1 mm∕h and above, however, the number density of scatterers is high enough to cause a deviation from the Beer-Lambert result. An increasing fraction of the radiation is interacting with the snow, and therefore an increasing amount of radiation will be scattered within the target solid angle, increasing the transmission efficiency compared to the Beer-Lambert result. It can be noted that the transmission efficiencies are not monotonically decreasing with increasing precipitation rate. As discussed below, the reason lies in variations in the PSDs and the classification of precipitation type. Figures 7 and 8 show comparisons between the experimental signal and the transmission efficiency calculations. The figures also show a curve fit of the experimental signal based on a typical power law assumption, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 2 7 2 η atm ∝ I g prec : In Eq. (11) η atm is the atmospheric transmission efficiency, I prec is the precipitation rate, and g is the fitting parameter. Figure 7 shows the experimental signal and the result based on transmission calculations using the Beer-Lambert transmission approach from Eq. (5), and Fig. 8 shows the experimental signal together with the result based on transmission calculations using the forward scattering transmission approach from Eq. (9). In addition, coefficients of determination, R 2 , indicating correlation between experimental and calculated data are presented in the figures as well. An R 2 -value equal to 100% corresponds to complete correlation between experiment and calculation, and an R 2 -value equal to 0% corresponds to the calculated signal being constant and equal to the mean value of the experimental signal. As can be seen in Figs. 7 and 8, both the Beer-Lambert approach and the approach including forward scattering perform significantly better in capturing the variation in measured signal than the curve fit based on the experimental signal and the recorded precipitation rate. Nonmonotonic variations in the signal with respect to the precipitation rate can be accounted for using the details of the PSDs as opposed to the precipitation rate only. Comparison of R 2 -values Fig. 6 Calculated transmission efficiency based on Mie theory and the disdrometer data as a function of snow precipitation rate. η at mBL , which is based only on the extinction coefficient result in the lowest transmission, whereas η at mF S−out give the highest transmission. Both transmission efficiencies based on inclusion of forward scattering result in higher transmission than the pure extinction case. The non-monotonous behavior of the calculated transmission efficiencies is caused by variations in the snow particle density and PSD.
of the approach using the Beer-Lambert law with the approach including the contribution from forward scattering show that inclusion of singly scattered light will provide better agreement with experiment than the Beer-Lambert approach. The comparatively large area of the target allow light scattered in the forward direction to contribute to the signal, increasing the outbound transmission efficiency compared to the Beer-Lambert approach, as shown in Fig. 6. The forward scattering approach only fails to capture the signal level at the highest precipitation rate of 2.44 mm∕h, whereas the Beer-Lambert approach start to underestimate the signal level already at a precipitation rate of 0.24 mm∕h. The reason for the underestimation of the signal level by the single scattering approach at 2.44 mm∕h cannot be conclusively determined since the precipitation characterization was performed with pointwise measurements and the transmission measurements along an extended path. In the discussion section, two hypotheses concerning spatially varying precipitation and multiple scattering, respectively, are discussed.
To exemplify the effects of PSD variation, two parts of the scaled signal versus precipitation rate curve have been investigated. The first part is the drop in signal at a fixed precipitation rate of   Fig. 9. It is evident that even though the two cases with lowest signal have significantly lower number of particles with diameter of 0.44 mm and smaller, the fact that they have larger number of particles with diameter of 0.6 mm and larger, make transmission lower here. Despite the lower number density, the larger diameter particles of the PSD affect transmission more than the small diameter particles. Interpretation of the difference between the cases with scaled signal values of 2.7 and 1.7 is straightforward. The case with a scaled signal value of 1.7 have a larger number of particles for all diameters except for the largest measured size, corresponding to a diameter of 2.75 mm, giving rise to a lower transmission. Despite the differences in PSD, the disdrometer has calculated an equal precipitation rate for the three cases. In the case of the scaled signal value of 5.8, the disdrometer classified the precipitation as a mix of snow grain and liquid drizzle. The other two cases, with scaled signal values of 2.7 and 1.7, were both classified as containing liquid drizzle, snow grains, and graupels, in different ratios. The classification of precipitation type is done by internal disdrometer algorithms based on temperature and correlation between velocity and size of the measured hydrometeor. The reason for the reported equal precipitation rates despite the difference in PSD is the mass density, which is different for drizzle, snow grains, and graupels.
The second part consists of the three adjacent cases corresponding to precipitation rates of 0.52, 0.67, and 0.90 mm∕h. The measured signal levels for these cases can be seen in Figs. 7 and 8, and they were chosen because of the clear reduction of both signal and calculated transmission in the 0.67 mm∕h case, when compared to the surrounding data points. The PSDs corresponding to the chosen cases are shown in Fig. 10. The factor responsible for the reduction in transmission efficiency is the higher number density in the diameter interval between 0.63 and 3.3 mm for the 0.67 mm∕h case. Within this range, the combination of snow particle number density and snow particle extinction cross-section generate a lower signal than in the other two cases, even though the precipitation rate placed this case in between the two others.

Discussion
Our results show that basing the transmission calculations on measured hydrometeor PSD and hydrometeor scattering behavior give a good agreement between calculations and measured signal levels in the case of falling snow. Furthermore, if the effective areas of the target and the detector are known, including forward scattering in the calculations will give a better estimation of the atmospheric transmission than using only the Beer-Lambert extinction approach. In our case, only at the highest recorded precipitation rate did the forward scattering approach fail to capture the signal level. Since the precipitation was measured at a single point whereas the transmission was measured along an extended path, we cannot conclusively determine the reason for this underestimation.
One hypothesis is that the underestimation of signal level at the highest precipitation rate is due spatially inhomogeneous precipitation. The precipitation was dominated by snow showers drifting by the measurement site, as can be inferred from the temporal variation shown in Fig. 4. It is thus likely that, to varying extent, there was spatially inhomogeneous precipitation for the captured measurement cases. This spatial variation is one hypothesis that can explain the difference between experimental transmission and calculations. Especially during the higher recorded precipitation rates, it is likely that the mean precipitation rate along the path is lower than the point measurement indicates. This would give rise to a higher transmission efficiency than calculated and can explain the underestimation of the calculated transmission efficiency at the highest recorded precipitation rate.
A second hypothesis that can explain the difference is effects from multiple scattering. The attenuation length (AL) is defined as B ext L, using the nomenclature from Eq. (5), and AL equal to 1 is often recognized as the onset of significant multiple scattering effects. The AL for the measured cases is shown in Fig. 11. Here it can be seen that AL is larger than 1 for precipitation rates of 0.24 mm∕h and higher. Thus, potentially, for this particular geometry and type of scattering particles, the single scattering approximation is valid for all but the highest recorded precipitation rate, where multiple scattering effects will start to have an impact.
Relying on empirical power laws with respect to the precipitation rate, measured with the laser disdrometer, when estimating the atmospheric transmission give significantly worse agreement with experimental data than using calculations based on the PSD. The large variation of transmission with respect to precipitation rate in falling snow have previously been noted in in the literature, and have been attributed to variations in snow crystal growth and aggregation, affecting the mass density and the PSD. 7,8 The correlation between the precipitation rate and both the measured transmission signal and the data from the visibility sensor, are in agreement with these assessments. However, as the work by Fehlmann et al. 15 showed, the operating principle of the laser disdrometer, relying on hydrometeor mass density estimations, will result in uncertainties when transforming the raw data to a precipitation rate in the case of falling snow. Therefore, it is difficult to assess if our observed large variations in transmission with respect to precipitation rate is due to inherent variations in the precipitation, or variations introduced by the measurement equipment, or both. Further measurements using a disdrometer for the PSD and a high accuracy precipitation rate measurement device, for example a weighing precipitation gauge, can give a more accurate understanding of the role of precipitation rates with respect to transmission of laser radiation in falling snow. The approximation of treating the snow PSD as composed of spherical objects is rather crude. 26,27 Even though, this approximation gives an improvement in agreement with experimental data compared to the extinction only approach. However, a better approximation of the scattering phase function will likely lead to a better agreement between calculations and experiments.
In the experimental setup, the target was much larger than the detector aperture. If instead a small target is used in a lidar geometry, or if a small detector is used in a single path setup, the difference between the forward scattering and extinction approaches will be smaller than what was measured in our experiments. The Beer-Lambert approach would then likely give good agreement with experiments for higher precipitation rates than in the investigated geometry.
The effect on imaging contrast and resolution due to forward scattering in the path from the target to the detector has not been included in the discussion here. Instead, all collected photons within the full image of the target have been summed. Photons scattered in a different GmAPD pixel than the one aimed at the part of the target where the photon was reflected will lead to image degradation. As seen from the small difference in η atmBL and η atmFS−in in Fig. 6, the number of photons scattered in the path from the target to the sensor that are collected by the detector is relatively small, and the image degradation should be limited. In high scattering environments, such as water with added scattering agents, taking an example from a different application, substantial effects have been measured. 28 Scattering in the outbound path from the laser to the target will not affect the resolution.

Summary
In this work, measurements of transmission in falling snow in a lidar geometry using a TCSPC 3D imaging lidar with a wavelength of 1.56 μm are presented. The snow PSD was characterized with a laser disdrometer providing a point-wise characterization of the precipitation, and the PSDs were used to calculate transmission efficiencies based on Mie scattering behavior.
The transmission was calculated in two ways, using the Beer-Lambert approach where all radiation interacting with a hydrometeor are considered extinct and an approach that includes effects of the forward scattering within the single scattering approximation. Comparison of these two methods showed that including the contribution from forward scattering gives better agreement between experiment and calculations than using extinction only. Including forward scattering effects allow signal behavior to be accurately captured in higher precipitation rates than the extinction approach. When comparing the results using details of the PSD with a curve fit approach based on precipitation rate, it is evident that both investigated approaches give a much better correlation between experiment and calculations than relying on precipitation rate only. The used disdrometer, however, has inherent uncertainties when converting raw data to precipitation rates. The reason is uncertainties in estimation of snow particle mass density. A higher accuracy precipitation rate gauge would shed further light into the role of precipitation rate with respect to transmission of laser radiation in falling snow.