Igor Veselovskii

Qiaoyun Hu
email: qiaoyun.hu@univ-lille.fr

Philippe Goloub

Thierry Podvin

William Boissiere

Mikhail Korenskiy

Nikita Kasianik

Sergey Khaykyn

Robin Miri

Derivation of depolarization ratios of aerosol fluorescence and water vapor Raman backscatters from lidar measurements

Polarization properties of the fluorescence induced by polarized laser radiation are widely considered in laboratory studies. In lidar observations, however, only the total backscattered power of fluorescence is analyzed. In this paper we present results obtained with a modified Mie-Raman-fluorescence lidar operated at the ATOLL observatory, Laboratoire d'Optique Atmosphérique, University of Lille, France, allowing us to measure depolarization ratios of fluorescence at 466 nm (δ F ) and of water vapor Raman backscatter. Measurements were performed in May-June 2023 during the Alberta forest fires season when smoke plumes were almost continuously transported over the Atlantic Ocean towards Europe. During the same period, smoke plumes from the same sources were also detected and analyzed in Moscow, at the General Physics Institute (GPI), with a five-channel fluorescence lidar able to measure fluorescence backscattering at 438, 472, 513, 560 and 614 nm. Results demonstrate that, inside the planetary boundary layer (PBL), the urban aerosol fluorescence is maximal at 438 nm, and then it gradually decreases with the increase in wavelength. The smoke layers observed within 4-6 km height present a maximum fluorescence at 513 nm, while in the upper troposphere, fluorescence maximum shifts to 560 nm. Regarding the fluorescence depolarization ratio, for smoke its value typically varies within the 45 %-55 % range.

The depolarization ratio of the water vapor Raman backscattering at 408 nm is shown to be quite low (2±0.5 %) in the absence of fluorescence because the narrowband interference filter (0.3 nm) in the water vapor channel selects only the strongest vibrational lines of the Raman spectrum.

As a result, the depolarization ratio at the water vapor Raman channel is sensitive to the presence of strongly depolarized fluorescence backscattering and can be used for the evaluation of the aerosol fluorescence contribution to measured water vapor mixing ratio.

Introduction

The possibility to measure the laser-induced fluorescence becomes an important added value to existing Mie-Raman lidars because fluorescence measurements provide new independent information about aerosol properties. Nowadays, the spectroscopic lidars based on 32-channel PMT combined with spectrograph proved the ability to measure the fluorescence spectrum [START_REF] Sugimoto | Fluorescence from atmospheric aerosols observed with a multichannel lidar spectrometer[END_REF][START_REF] Reichardt | Cloud and aerosol spectroscopy with Raman lidar[END_REF]Reichardt et al., 2018Reichardt et al., , 2023;;[START_REF] Richardson | Bioaerosol detection over Athens, Greece using the laser induced fluorescence technique[END_REF][START_REF] Liu | Spectrally Resolved Raman Lidar to Measure Backscatter Spectra of Atmospheric Three-Phase Water and Fluorescent Aerosols Simultaneously: Instrument, Methodology, and Preliminary Results[END_REF]. On the other hand, lidars, with a single fluorescence channel, can be widespread due to their simplicity [START_REF] Rao | Preliminary measurements of fluorescent aerosol number concentrations using a laser-induced fluorescence lidar[END_REF]Veselovskii et al., 2020). Such singlechannel fluorescence lidars, combined with depolarization measurements at the elastic wavelength, provide new independent information about aerosol type (Veselovskii et al., 2022). However, in all lidar studies, only the total scattered power was analyzed, while the polarization properties of the fluorescence were ignored. At the same time, fluorescence depolarization measurements are widely used in laboratory research [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. When polarized laser radiation is used for excitation, the fluorescence emission is also partly polarized and the degree of its depolarization Published by Copernicus Publications on behalf of the European Geosciences Union.

(anisotropy) depends on the fluorescence lifetime, on the angle between excitation and emission dipoles, and on the rotational mobility of molecules [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. In the fluorescence spectroscopy, the polarization state of emission is described by the anisotropy [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF], introduced as

r = P F -P ⊥ F P F + 2P ⊥ F , (1)
where P F and P ⊥ F are the powers of co-and cross-polarized fluorescence components. In lidar measurements, however, the fluorescence depolarization ratio, δ F , is given as

δ F = P ⊥ F P F . (2)
Therefore, the anisotropy is expressed as a function of the δ F as follows:

r = 1 -δ F 1 + 2δ F . (3)
For randomly oriented fluorophores with collinear absorption and emission dipoles, in the absence of rotational motion, the anisotropy r = 0.4 [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF], which corresponds to δ F = 33 %. This is the minimal value one can expect in lidar measurements. The existence of any angle between absorption and emission dipoles, as well as molecule rotation in the process of emission, will increase δ F [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. Thus, measurement of fluorescence depolarization ratio may bring additional information about atmospheric aerosol, as we will show below.

Water vapor is a key atmospheric component playing an essential role in the planet's radiative balance, and Raman lidars today are widely used for such observations [START_REF] Whiteman | Examination of the traditional Raman lidar technique. I. Evaluating the temperature dependent lidar equations[END_REF]Chouza et al., 2022, and references therein). However, when the UV laser beam passes through a smoke layer, the broadband fluorescence signal is induced and its spectrum includes the region of water vapor Raman lines. Thus, the signal in the water vapor channel (around 407.5 nm, when 354.7 nm laser radiation is emitted) becomes contaminated by the fluorescence backscatter signal [START_REF] Immler | Is fluorescence of biogenic aerosols an issue for Raman lidar measurements?[END_REF][START_REF] Immler | Is fluorescence of biogenic aerosols an issue for Raman lidar measurements?[END_REF]. This contamination can be reduced by decreasing the width of the transmission band in the water vapor channel down to 10ths of nanometers. However, as was shown recently, fluorescence still remains the source of uncertainties, especially when the water vapor mixing ratio (WVMR) is measured inside the smoke layers in the upper troposphere (Chouza et al., 2022;Reichardt et al., 2023).

Depolarization measurements provide an opportunity to monitor the presence of fluorescence signals in the Raman channel. The Q branch of water vapor Raman lines (near 407.5 nm) provides a weakly depolarized backscatter, while fluorescence is strongly depolarized. Thus, the presence of fluorescence should increase the depolarization ratio of the signal in the water vapor channel. Moreover, if the depolarization ratios of water vapor and fluorescence are known, the contribution of fluorescence to the measured WVMR can be evaluated.

In this article, we report and analyze, for the first time, the depolarization ratio of aerosol fluorescence and of water vapor Raman backscatter from lidar observations performed at the ATOLL observatory (ATmospheric Observation at liLLe), Laboratoire d'Optique Atmosphérique, University of Lille, during dense smoke events occurring in May-June 2023. We start with a description of the experimental setup in Sect. 2.1 and derive, in Sect. 2.2, the main equations for estimating the fluorescence contribution to the water vapor Raman channel. In the first part of the results section (Sect. 3.1), the fluorescence depolarization ratios over ATOLL are analyzed for different aerosol types. The measurements of fluorescence spectra performed with a new five-channel fluorescence lidar, operated in Moscow, are presented in Sect. 3.2. In Sect. 3.3, we analyze the depolarization ratio in the water vapor Raman channel and estimate the contamination of fluorescence in the derived WVMR profiles. Finally, in Sect. 4 we present our conclusions.

2 Experimental setup and data analysis

Lidar system

In our study, two lidar systems are considered. The first one, LILAS (LIlle Lidar AtmosphereS), is a multiwavelength Mie-Raman-fluorescence lidar, whereas the second one is a multiwavelength fluorescence lidar operated by the General Physics Institute (GPI), Moscow (Veselovskii et al., 2023). Both systems are based on a tripled Nd:YAG laser (Q-smart 450) with a 20 Hz repetition rate and pulse energy about of 100 mJ at 355 nm. The backscattered laser light in both systems is collected by a 40 cm aperture telescope, and the lidar signals are digitized with transient recorders (Licel) with 7.5 m range resolution, allowing simultaneous detection in the analog and photon counting modes.

LILAS allows the so-called 3β + 2α + 3δ configuration, including three particle backscattering (β 355 , β 532 , β 1064 ), two extinction (α 355 , α 532 ) coefficients, and three particle depolarization ratios (δ 355 , δ 532 , δ 1064 ). The Raman channel with a 407.54/0.3 nm spectral width interference filter allows also water vapor profiling. At the end of 2019, the lidar was modified to enable fluorescence measurements. A part of the fluorescence spectrum is selected by a wideband interference filter of 44 nm width centered at 466 nm (Veselovskii et al., 2020).

In the fluorescence lidar of GPI only a 355 nm wavelength is emitted, while fluorescence is measured in five spectral intervals. The central wavelengths and widths of spectral transmission bands (in parentheses) are 438 (29), 472 (32), 513 (29), 560 (40) and 614 (54) nm (Veselovskii et al., 2023). Thus, the fluorescence spectrum could be sampled at five different wavelengths. The transmission bands of the fluorescence channels (Fig. 1 in Veselovskii et al., 2023) are separated, and there are no cross-talks between the channels. At GPI, the measurements were performed at an angle of 48 • to the horizon. The strong sunlight background restricts the fluorescence observations of both systems to only the nighttime hours.

Several aerosol properties can be derived from fluorescence. The fluorescence backscattering coefficient, β Fλ , at wavelength λ F , is calculated from the ratio of fluorescence and nitrogen Raman backscattering signals, as described in Veselovskii et al. (2020). Remember that β Fλ is related to fluorescence signals integrated over the filter transmission band D λ . In Moscow measurements are performed at five wavelengths, and to compare β Fλ between different channels one makes use of the "fluorescence spectral backscattering coefficient" B λ = β Fλ D λ (fluorescence backscattering per spectral interval). LILAS has only one single fluorescence channel; therefore, when presenting data from LILAS, for the sake of simplicity, one uses the notation β F 466 = β F . The intensive property characterizing aerosol fluorescence is the fluorescence capacity G Fλ , which is the ratio of the fluorescence backscattering at wavelength λ F to the backscattering coefficient at laser wavelength G Fλ = β F β λ . This ratio, in principle, can be calculated for any laser wavelength. For LILAS observations G Fλ is calculated with respect to β 532 , as β 532 is derived with rotational Raman scattering and it does not depend on assumption about the Ångström exponent (Veselovskii et al., 2015). Again, when presenting LILAS data, for simplicity one will use the notation G Fλ = G F . In this work, all profiles of aerosol properties are smoothed with the Savitzky-Golay method, using second-order polynomials with eight points in the spatial window.

Additional information about the atmospheric thermodynamic state was available from radiosonde measurements performed at Herstmonceux (UK) and Beauvechain (Belgium) stations, located 160 and 80 km away from the ATOLL observatory, respectively. When calculating the relative humidity, the water vapor profiles measured by Raman lidar and temperature profiles provided by the radiosonde were used.

As discussed in Sect. 1, measurements of the fluorescence depolarization ratio and the depolarization of water vapor Raman backscatter are expected to bring new information about aerosol properties and fluorescence contamination in the water vapor Raman channel. In 2023, LILAS was upgraded to allow depolarization measurements at both 466 and 408 nm. The corresponding optical layout is shown in Fig. 1. Dichroic mirrors (DMs) separate the 387, 408 and 466 nm components, while polarizing cubes split the components with polarizations oriented parallel and perpendicular to the emitted polarized laser beam. For both channels, the polarizing cube PBS251 from Thorlabs was used. The fluorescence depolarization ratio, δ F , and the water vapor Raman scattering depolarization ratio, δw, are both defined and cal- culated as a ratio of the perpendicular to the respective parallel components. The calibration of both ratios was performed as described in [START_REF] Freudenthaler | Depolarization ratio profiling at severalwavelengths in pure Saharan dust during SAMUM[END_REF]. The uncertainty of calibration is estimated to be below 15 % for both 466 and 408 nm channels.

Expressions for estimating fluorescence impact on water vapor measurements

As discussed in the recent work of Chouza et al. (2022) and Reichardt et al. (2023), the broadband aerosol fluorescence is expected to contribute to the signal measured at the water vapor Raman channel. Below, we provide the basic equations for estimating this contribution, based on the measurements of the depolarization ratio in the water vapor Raman channel. The elastic backscattered radiative power, at the laser wavelength λ L from distance z, can be modeled, after background subtraction, by writing the following lidar equation:

P L = O(z) 1 z 2 C L βT 2 L , (4)
where O(z) is the geometrical overlap factor, which is assumed to be the same for all channels. C L is a rangeindependent constant, including efficiency of the detection channel, the emitted laser power and the receiving telescope diameter. T L is the one-way atmospheric transmission, describing light losses on the way from the lidar to distance z at wavelength λ L .

T L = exp    - z 0 α a (λ L , z ) + α m (λ L , z ) dz    (5)
The backscattering and extinction coefficients contain the aerosol (a) and molecular (m) contributions:

β λ L = β a λ L +β m λ L and α λ L = α a λ L + α m λ L .
Radiative power in nitrogen Raman, water vapor Raman and fluorescence channels can be written in a similar way.

P R = O(z) 1 z 2 C R σ R N R T L T R , (6)
https://doi.org/10.5194/amt-17-1023-2024 Atmos. Meas. Tech., 17, 1023-1036, 2024

P W = O(z) 1 z 2 C W N W σ W T W T L , (7)
P F = O(z) 1 z 2 C F β F T F T L , (8)
where C R , C W and C F are the corresponding rangeindependent constants. T R , T W and T F are the one-way transmissions at wavelengths λ R , λ W and λ F corresponding to the centers of transmission bands of the channels. N R and N W are the concentrations in nitrogen and water vapor molecules, while σ R and σ W are their Raman differential scattering cross sections, respectively. The fluorescence backscattering coefficient, β F , is introduced the same way, as described in Veselovskii et al. (2020).

The received power of the fluorescence signal that leaks to the water vapor channel is

P FW = O(z) 1 z 2 C W β FW T W T L , (9)
where β FW is the fluorescence backscattering coefficient at wavelength λ W . The WVMR, n W , can be obtained from Eqs. ( 6) and ( 7) if the calibration constant

K W = C R C W σ R σ W is known: n W = K W P W P R T R T W . (10)
The fluorescence backscattering coefficient, β F , derived from Eqs. ( 6) and ( 8), also contains the calibration constant K F . The procedure of calibration is described in Veselovskii et al. (2020). Finally, β F reads as

β F = K F n R P F P R T R T F , (11)
where z=0) is the relative change in number density of nitrogen molecules with height.

n R = N R (z) N R (
The fluorescence signal P FW in the water vapor channel can be expressed from P F using parameter η, which depends on the ratio of fluorescence cross sections at wavelengths λ W and λ F , on the filter's width, and on the efficiency of both channels, as follows:

P FW = P F η T W T F . ( 12
)
The total signal measured in the water vapor channel, PW , is the addition of both water vapor backscatter, P W , and the fluorescence backscatter, P FW :

PW = P W + P FW = P W + P F η T W T F . ( 13
)
One should remember that the fluorescence spectrum, even for the same type of aerosols, can vary with altitude and from observation to observation, which finally influences η.

To minimize this influence it is desirable to keep λ W and λ F as close as possible.

If the received lidar signals at the water vapor Raman and fluorescence channels are separated into co-polarized ( ) and cross-polarized (⊥) components, with respect to the polarization of the emitted laser beam, their powers at the water vapor Raman channel are given, respectively, as

P W = P W + P F η T W T F , (14)
P ⊥ W = P ⊥ W + P ⊥ F η T W T F = δ W P W + δ F P F η T W T F , (15)
where δ F and δ W are the fluorescence and water vapor Raman depolarization ratios, defined as

δ F = P ⊥ F P F and δ W = P ⊥ W P W . ( 16
)
Here we assume that the depolarization ratio of fluorescence is the same at the wavelengths λ W and λ F . This assumption is usually valid because fluorescence emission is normally from the lowest singlet state, so the depolarization ratio is spectrally independent [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. Due to the presence of fluorescence, the depolarization ratio measured at the water vapor Raman channel is

δW = P ⊥ W P W = δ W P W + δ F P F η T W T F P W + P F η T W T F . ( 17
)
Here δ W is the depolarization ratio that would be measured at the water vapor Raman channel in the absence of atmospheric fluorescence. From Eqs. ( 9), ( 10), ( 14), ( 15) and ( 17) the parameter η can be derived using the lidar-measured values, such as the water vapor mixing ratio ñW , depolarization ratio δW and fluorescence backscattering β F :

η = ñW β F K F K W n R (1 + δ F ) δW -δ W 1 + δW (δ F -δ W ) , ( 18
)
where ñW is the WVMR containing the fluorescence contribution.

It should be noted that the choice of calibration constants K F and K W does not influence η because ñW and β F are calculated using the same calibration constants. Finally, the increase in WVMR n W induced by the fluorescence can be calculated as

n W = K W P F η T W T F P R T R T W = K W K F ηβ F 1 n R . ( 19
)
As soon as the parameter η is calculated from Eq. ( 18), we can estimate the relevant error n W from β F , which in the case of LILAS is measured at 466 nm (Veselovskii et al., 2020). In such estimation we have to assume that the relationship between fluorescence at 466 and 408 nm remains constant with height. A possibility to perform correction from single-channel fluorescence measurements was discussed by Reichardt et al. (2023), where it was shown that for 466 and 408 nm channels, the correction actually may depend on height. The corresponding analysis based on our measurements will be presented in Sect. 3.2. We should mention that when the depolarization at the water vapor Raman channel is available, the contribution of fluorescence to WVMR can be obtained without using η. From Eqs. ( 18) and ( 19) we obtain

n W = ñW (1 + δ F ) δW -δ W 1 + δW (δ F -δ W ) . ( 20
)
However, such correction can be performed only at low altitudes, where the signal-to-noise ratio at the cross-polarized water vapor channel is sufficient for the calculation of δW .

Experimental results

In May-June 2023, the Canadian forest fires were the origin of numerous smoke layer observations in a wide range of altitudes, ranging from the planetary boundary layer (PBL) to the tropopause. The boreal wildfire season in 2023 started anomalously early. A wildfire in Alberta, Canada, at 53.2 • N, 115.7 • W produced an intense pyrocumulonimbus (PyroCb) cloud on 5 May with the minimum satellite-derived infrared brightness temperature of -66 • C, which should correspond to 10-11 km altitude according to local radiosoundings. In order to describe the long-range transport of the smoke plume produced by this event, we use UV absorbing Aerosol Index (AI) measurements by the Ozone Monitoring and Profiling Suite (OMPS) Nadir Mapper (NM) instrument on board the Suomi National Polar-orbiting Partnership (NPP) satellite mission (Flynn et al., 2014). AI is widely used as a proxy for the amount of absorbing aerosols (e.g., smoke, dust, ash), and its dimensionless value is proportional to the altitude of the aerosol layers. AI values above 15 are usually associated with smoke plumes at or above the tropopause (Peterson et al., 2018, and references therein), whereas the maximum AI value reported by the OMPS NM instrument for the Alberta event reached a value of 19.9. Figure 2 displays the spatiotemporal evolution of the smoke plume from the Alberta event represented by the areas of enhanced AI observed between 5 and 21 May. The smoke in the upper troposphere and lower stratosphere (UTLS) is advected by the westerly winds, crossing the Atlantic about 1 week before reaching Moscow on 15 May. On that date, the Moscow lidar detected the smoke layer at 10-11 km (see Sect. 3.2). The plume was then further advected across Eurasia towards northeastern Siberia. By 22 May the smoke plume completed its first circumnavigation (not shown) and passed over Lille on 23 May and then over Moscow for the second time around 27 May. Thus, we can expect that the smoke layers observed over Lille and Moscow have the same source.

Variability in fluorescence depolarization ratio

At the first stage of our research we focused on the variability in the fluorescence depolarization ratio for aerosol types. The main attention was paid to smoke particles because they provide the strongest impact on the water vapor Raman measurements due to their high fluorescence capacity (Veselovskii et al., 2022).

Spatiotemporal distributions of the aerosol elastic and fluorescence backscattering coefficients (β 532 and β F ), on the night of 26-27 June 2023, are shown in Fig. 3. A dense smoke layer with β F as high as 7.0 × 10 -4 Mm -1 sr -1 occurred within the 4.0-10.0 km height range. The HYSPLIT back trajectories show that the air masses were transported from North America. The relative humidity increased from 40 % at 4 km to > 90 % at 7 km, where the formation of ice crystals started. Vertical profiles of aerosol elastic and fluorescence backscattering coefficients (β 532 and β F ), together with fluorescence capacity, are shown in Fig. 3c. Inside the smoke layer, G F is about 3 × 10 -4 , which is a typical value for smoke (Veselovskii et al., 2022), whereas, above 6 km, it decreases due to ice formation. The presence of ice crystals increased the particle depolarization ratio δ 532 from 3 % at 6 km to 20 % at 8 km. Fluorescence signals are strongly depolarized. Inside the PBL, δ F was about 60 % whereas above 2 km it dropped to approximately 45 %. The processes of hygroscopic growth and ice formation do not provide a noticeable impact on the δ F value. During May-June observations, the depolarization ratio of smoke varied mainly inside the 45 %-55 % range.

As discussed in our previous publications (Veselovskii et al., 2022;Hu et al., 2022), the fluorescence capacity of aged smoke varies inside the (2.5-5.5) ×10 -4 range, probably due to the changes in smoke composition and conditions of atmospheric transport. However, during the Alberta fires, several smoke plumes with high G F have been observed. The highest fluorescence capacity was observed on the night of 16-17 June 2023. Vertical profiles of the aerosol properties for this episode are shown in Fig. 4. A dense smoke layers with fluorescence backscattering exceeding 10.0 × 10 -4 Mm -1 sr -1 occurred within the 7.0-9.0 km height range. In this case, the maximal value of the fluorescence capacity reached 10.0 × 10 -4 . The fluorescence depolarization ratio was measured about 50 % through the entire smoke layer, and the process of ice formation (just like in Fig. 3d) does not influence δ F . Thus, in May-June 2023 strong variations in G F in the (2.5-10.0) ×10 -4 range were accompanied by relatively small variations in δ F remaining in the 45 %-55 % interval.

It is known that in the UTLS smoke particles can reach a depolarization ratio, δ 532 , as high as 15 %-20 % ( , 2020). High values of the particle depolarization ratio are usually attributed to the complex internal structure of smoke particles [START_REF] Mishchenko | Linear depolarization of lidar returns by aged smoke particles[END_REF]. Two smoke events in the UTLS, characterized by enhanced δ 532 on 28-29 May and 3-4 June 2023, are illustrated in Fig. 5. On 28-29 May, three smoke layers, at ∼ 3.5, 6.5 and 11.5 km, can be distinguished. High depolarization ratios, reaching 40 % at altitudes of 9.8-10.5 km, are due to ice clouds. In the lower smoke plumes ranging between 3.5 and 6.5 km, the particle depolarization did not exceed 8 %, whereas above 11 km δ 532 increased to 15 %. High values of δ 532 observed in the UTLS correlate with an increase in G F and with fluorescence depolarization, δ F , up to 7.0 × 10 -4 and 70 %, respectively. Similar behavior was observed on 3-4 June, when the depolarization ratio, δ 532 , above 11.5 km increased up to 15 %, simultaneously with an increase in G F and δ F up to 9.5 × 10 -4 and 70 %, respectively. Thus, the change in particle morphology may affect the depolarization ratio at the fluorescence channel. Another possibility is that, in the UTLS, not only the particle structure can change, but the composition can as well. At the current stage of analysis, we are not yet able to draw conclusions about the mechanisms explaining the increase in fluorescence depolarization in the UTLS. Furthermore, we did not observe the effect of atmospheric humidity on smoke fluorescence depolarization. However, inside the PBL the observed hygroscopic growth was accompanied by an increase in δ F . During the 9-16 June 2023 period numerous particle hygroscopic growth cases were ob- served in the PBL. One such case, on the night of 12-13 June, is shown in Fig. 6. The relative humidity increased inside the PBL from 50 % to 90 % causing an increase in β 532 near the PBL top. Depolarization ratio δ 532 decreased with height, since the particles in the process of hygroscopic growth became more spherical. The fluorescence depolarization ratio, however, increased inside the PBL from 50 % to 70 %.

All results obtained during 9-16 June, showing the dependence of δ F and δ 532 on the relative humidity, are summarized in Fig. 7. Particle depolarization δ 532 systematically decreased with relative humidity (RH), but, on 16 June, this dependence was not monotonic, which could be due to the change in aerosol composition with height. At low RH (below 30 %), the fluorescence depolarization ratio was about 50 %. However, at RH of about 90 %, δ F increased up to 70 %. One of the possible explanations for that behavior could be an increase in rotational mobility of the molecules in the process of particle water uptake.

Fluorescence spectrum sampled with a five-channel lidar

The results presented in the previous section were obtained with a single-channel fluorescence lidar. However, for analyzing the variability in smoke properties (for example, increase in the fluorescence capacity with height), it is important to have information about a wider fluorescence spectrum. Moreover, to estimate the fluorescence contamination in the water vapor Raman channel, a relationship between fluorescence backscattering at 466 and 408 nm is used. Thus, we need to know the variability in the fluorescence spectrum in the short wavelength region. In our recent work (Veselovskii et al., 2023) we presented the first results obtained with a five-channel fluorescence lidar in operation at the GPI. This lidar is able to measure the fluorescence backscattering profiles at five spectral intervals centered at 438, 472, 513, 560 and 614 nm. During May-June 2023, several smoke plumes originating from Alberta fires were transported over Moscow. Although Lille and Moscow are very distant from each other (above 2200 km), the smoke plumes observed have the same origin; hence the fluorescence spectra measured over Moscow are quite helpful for the analysis of the Lille data. Figure 8a, b and c present the fluorescence spectral backscattering coefficients, B λ , for three smoke events detected in the UTLS above 10, 8 and 10 km for 15 May, 31 May and 20 June 2023, respectively. On 15 and 31 May smoke layers were also present inside the 4-6 km range. Inside the PBL the strongest fluorescence was systematically detected at the 438 nm channel, while, at higher altitudes, the maxima shifted to 560 nm. As follows from Figs. 8df, the ratio B 560 /B 438 remained in the range 0.4-0.7 inside the PBL, whereas this ratio increased above 2.0 in the UTLS. Thus, for smoke events the maxima of the fluorescence spectrum shifted with height towards longer wavelengths. The ratio B 513 /β 355 also increased with height, and, above 10 km, it reached the values of 1 × 10 -5 nm -1 . In the UTLS, the maximal fluorescence capacity, G F , measured by LILAS at 466 nm (with 44 nm bandwidth filter) was about 10 × 10 -4 . In the smoke layer, the ratio of backscattering coefficients β 355 /β 532 is about 2, so the maximal ratio B 466 /β 355 derived from LILAS measurements was about 1.1 × 10 -5 nm -1 . Thus, values obtained over Lille and over Moscow are in good agreement.

The fluorescence spectra obtained for the abovementioned smoke plumes are shown in Fig. 9. The values of B λ are normalized to B 438 . Inside the PBL, the maximum of fluorescence was measured at 438 nm, and it decreased with wavelength. In the smoke layers within 4-6 km, the maximum of fluorescence is observed at 513 nm, while, in the UTLS, the maximum shifted to 560 nm.

When applying Eq. ( 19) to estimate the contribution of smoke fluorescence into the water vapor Raman channel of LILAS, we assumed that the ratio of the fluorescence backscattering at 466-408 nm (B 466 /B 408 ) was constant. For the lidar in operation at GPI, the shortest available wavelength was 438 nm; therefore, at least, one can estimate the variability in the ratio B 472 /B 438 . Figure 10 presents the vertical profiles of B 472 /B 438 for 11 smoke events occurring during the 15 May-20 June 2023 period. Inside the PBL, this ratio varied in the 0.6-1.0 range. The lowest values correspond to urban aerosols, while values of B 472 /B 438 close to 1.0 probably indicate the presence of smoke particles inside the PBL. Smoke layers were observed mainly above 4.0 km and B 472 /B 438 showed a tendency to increase in the UT. It is interesting that, for the period 15 May-1 June, the ratio was close to 1.5, whereas after 1 June, it became close to 1.0, which can be related to changes in the smoke source. The mean value of B 472 /B 438 in the 4.0-11.0 km range over all https://doi.org/10.5194/amt-17-1023-2024 Atmos. Meas. Tech., 17, 1023-1036, 2024 observations is 1.38 with a standard deviation of 0.23 (relative variation is about 17 %). For the 466 and 408 nm channels the wavelength separation is larger, so one can expect a variation in B 466 /B 408 in the smoke layer to be above that value. It points out the difficulties faced when the estimation of the fluorescence contamination in the water vapor Ra-man channel is performed from a single fluorescence channel at 466 nm. This issue was also discussed by Reichardt et al. (2023).

Estimation of fluorescence impact on water vapor Raman measurements

Measuring the depolarization ratio at the water vapor Raman channel provides an opportunity to control/evaluate the presence of a fluorescence leak in this channel. These depolarization measurements were performed over Lille during May-June 2023. Vertical profiles of water vapor depolarization ratio δW , together with ñW , β 532 , β F and G F , are shown in Fig. 11 for the nights of 8-9 and 10-11 June 2023. On 8-9 June the aerosols were confined mainly below 5 km. The fluorescence capacity was about 1.0×10 -4 below 3.0 km, but above, G F increased up to 2.5×10 -4 , indicating the presence of smoke particles. The depolarization ratio in the water vapor channel was about 2 % in the height range 1.5-3.5 km, where the values of δW ranging within 1.8 %-2.0 % were observed at this height range, where the contribution of fluorescence was insignificant. The depolarization ratio δ W was low because the interference filter at the water vapor channel selects only the strongest Q-branch lines and most of the rotational lines are blocked. The contribution of fluorescence becomes noticeable above 3.5 km where n W dropped, resulting in an increase in δW of up to ∼ 3 %. Below 1 km height we also observed an increase in δW up to 2.2 %, where fluorescence backscattering is enhanced. Similar values of δW were observed on 10-11 June, where the depolarization ratio increased up to 2.5 % inside the smoke layer observed at ∼ 3.75 km and below 2.0 km. As discussed in Sect. 2.2, the contribution of fluorescence to the WVMR channel can be derived from Eq. ( 20) if δW and δ F are measured simultaneously. Figure 12 presents the modeling of the relative error n W ñW introduced by the fluorescence to the WVMR channel as a function of δW . The computations are performed for different fluorescence depolarization ratios δ F = 50 %, 60 % and 70 % to include both smoke and urban particles. A depolarization ratio in the water vapor Raman channel in the absence of fluorescence was assumed to be δ W = 2 %. For a depolarization ratio δW below 3 % the relative error n W ñW did not exceed 3 %. As follows from the fluorescence spectra in Fig. 9, the fluorescence of urban particles increases towards shorter wavelengths; thus, one can expect an impact of the urban aerosol fluorescence on the water vapor measurements. In practice, however, we did not observe values of δW exceeding 3 % in the PBL; thus, the contribution of aerosol in the PBL is not critical. The reason is due to the low fluorescence capacity (about one order lower than that of smoke) and higher water vapor content compared to the free troposphere.

Vertical profiles of δW shown in Fig. 11 become noisy at heights where n W is low, and thus δW cannot be used for the correction of the fluorescence effect in the upper troposphere. To overcome this, we derived the parameter η from Eq. ( 18) at low altitudes where δW values are available, and, thus, these η values can be used to calculate n W from Eq. ( 19) in the entire height range. In such an approach, however, one has to assume that the relationship between fluorescence cross sections at 466 and 408 nm remains constant with height. As discussed in Sect. 3.2, such an assumption can yield significant bias in the calculation of n W , and, at this stage, we do not provide corrected profiles of the WVMR.

For the accurate calculation of η one needs smoke events with strongly enhanced δW values, which are usually observed in the smoke layers with low WVMR. Such suitable events are shown on the nights of 26-27 May and 5-6 June 2023 in Fig. 13. On 26-27 May a smoke layer characterized by high fluorescence (β F up to 5×10 -4 Mm -1 sr -1 ) and low ñW (below 0.2 g kg -1 ) values is observed at 3.5 km. The relevant fluorescence depolarization ratio was about 47 %, and δW increased from 2 % up to 12 % in the middle of this layer. The parameter η calculated from Eq. ( 18) inside this smoke layer was about 2×10 -3 (g kg -1 ) (Mm -1 sr -1 ) -1 . On 5-6 June the depolarization ratio δW in the smoke layer increased up to 10 %, and the value of η was very similar. The values of η derived for several smoke episodes varied in the range (2-2.5) ×10 -3 (g kg -1 ) (Mm -1 sr -1 ) -1 . To estimate n W we used the mean value of η = 2.25 × 10 -3 (g kg -1 ) (Mm -1 sr -1 ) -1 , which is suitable only for smoke, while for particles in the PBL, η can have a different value. However, in the PBL, the low depolarization ratios of δW prevented us from calculating η.

Figure 14 presents the vertical profiles of WVMR, the fluorescence backscattering and the error n W introduced by the fluorescence in WVMR on 26-27 May, 28-29 May and 16-17 June. Smoke layers with strong fluorescence occurred systematically in our upper-tropospheric observations. The current LILAS system is not powerful enough to derive accurate water vapor measurements above 10 km; however, an increase in ñW in the fluorescent smoke layers is visible. Remember that Eq. ( 19) for n W contains the factor 1 n R (inverse relative change in nitrogen number density); thus, the fluorescence impact on WVMR will increase with height. The uncertainties n W ñW for all events considered are shown in Fig. 14d. On 26-27 and 28-29 May the uncertainty of n W ñW at 11 km is of the order of 100 %. On 16 June the smoke layer is lower (at 9 km), and the uncertainty is about 50 %. Our demonstration shows that smoke fluorescence can significantly impact the water vapor measurements. The proposed approach, based on the analysis of the depolarization ratio of the water vapor signal, has the potential for the estimation and correction of this impact.

Conclusion

This study is one of the first efforts to measure the depolarization ratio of the fluorescence of the atmospheric aerosols. Analysis of more than 30 spring and summer smoke events allows the evaluation of the main aerosol intensive properties, including fluorescence capacity and the particle and fluorescence depolarization ratios. The fluorescence capacity of smoke in the troposphere varied within (2.5-10.0) ×10 -4 ; however, in spite of strong G F variation, δ F remained within a relatively narrow interval of 45 %-55 %. Additional observations revealed that for smoke plumes in the upper troposphere the fluorescence depolarization ratio increased up to 70 %. At the moment, we cannot fully explain the mechanism responsible for this δ F increase. It can be related to complex particle internal structure at high altitudes, as well as to the change in the chemical composition, revealed by the shift in the maximum of the fluorescence spectra to longer wavelengths in the upper troposphere (Fig. 9).

Inside the PBL, the fluorescence depolarization ratio was higher than that of smoke and varied within the 50 %-70 % range. Moreover, the fluorescence depolarization ratio of urban particles strongly depends on the relative humidity, and, in contrast to the elastic scattering, the depolarization of fluorescence increases with RH. One possible origin of this phenomena could be attributed to an increase in the rotational mobility of the molecules involved in the process of water uptake.

The depolarization ratio of the water vapor Raman backscatter, in the absence of fluorescence, appeared to be quite low (δ W = 2 ± 0.5 %). As a result, the depolarization ratio of the water vapor Raman backscatter is sensitive to the presence of strongly depolarized fluorescence signals, and the contribution of fluorescence to the WVMR can be calculated from the measured value δW . However, with the lidar used in this work, measurements of δW are only possible up to the middle troposphere, while the problem of the fluorescence interference is the most crucial in UTLS. To estimate the impact of fluorescence on the WVMR in UTLS, the height-independent parameter η, linking fluorescence at 466 and at 408 nm, was used. Such an approach relies on the assumption that η remains constant and allows only a rough estimation of the correction term for the WVMR, n W . One possible solution to increase the accuracy of n W is to implement an additional shorter wavelength channel (438 nm or even shorter). Another technical approach worth considering is equipping the 408 nm channel with a polarizing cube, given the low depolarization ratio of Raman water vapor backscatter. Thus, the depolarized channel at 408 nm can be used for fluorescence measurements. As the polarizing cubes work in a wide spectral range, one can select a spectral region outside of the water vapor spectrum (400-418 nm) for fluorescence monitoring. We plan this experiment, as well as other innovative approaches, with our future high-power fluorescence lidar, LIFE (Laser Induced Fluorescence Explorer), whose start of operation is scheduled for the beginning of 2024.

Figure 1 .

 1 Figure 1. Optical layout of depolarization measurements at 408 nm and 466 nm wavelengths. L -lens; IF1-IF3 -interference filters; DM1 and DM2 -dichroic mirrors; PBS -polarizing cube.

Figure 2 .

 2 Figure 2. Spatiotemporal evolution of the smoke plume from the wildfire event in Alberta, Canada, on 5 May 2023. Color-filled time-coded areas indicate the Aerosol Index (AI) values from the OMPS NPP instrument exceeding 0.5. The blue-and red-filled stars indicate the location of Lille and Moscow lidar stations, respectively.

Figure 3 .

 3 Figure 3. Smoke event on the night of 26-27 June 2023 over Lille. Spatiotemporal distributions of (a) aerosol backscattering coefficient β 532 and (b) fluorescence backscattering β F (in 10 -4 Mm -1 sr -1 ). Vertical profiles of (c) the aerosol β 532 and fluorescence β F backscattering coefficients and the fluorescence capacity G F , as well as (d) the particle δ 532 and the fluorescence δ F depolarization ratios.

Figure 4 .

 4 Figure 4. Vertical profiles of (a) aerosol β 532 and fluorescence β F backscattering coefficients, fluorescence capacity G F , and (b) particle δ 532 and fluorescence δ F depolarization ratios on the night of 16-17 June 2023 for the period 01:00-02:30 UTC over Lille.

Figure 5 .

 5 Figure 5. Vertical profiles of (a, c) backscattering coefficient β 532 , particle depolarization ratio δ 532 and fluorescence capacity G F , as well as (b, d) fluorescence backscattering β F and fluorescence depolarization ratio δ F , for two smoke episodes on the nights of 28-29 May and 3-4 June 2023 over Lille.

Figure 6 .

 6 Figure 6. Particle hygroscopic growth in the PBL on the night of 12-13 June 2023 over Lille. Spatiotemporal distributions of (a) aerosol backscattering coefficient β 532 and (b) fluorescence backscattering β F (in 10 -4 Mm -1 sr -1 ). Vertical profiles of (c) aerosol β 532 and fluorescence β F backscattering coefficients, (d) particle depolarization ratio δ 532 and the relative humidity RH, and (e) fluorescence depolarization ratio δ F and fluorescence capacity G F for the time period 21:00-22:15 UTC.

Figure 7 .

 7 Figure 7. (a) Fluorescence depolarization ratio and (b) particle depolarization ratio δ 532 as a function of the relative humidity in the PBL for the measurements on 9, 10, 12, 13 and 16 June 2023 over Lille.

Figure 8 .Figure 9 .

 89 Figure 8. Fluorescence measurements over Moscow on 15 May, 31 May and 20 June 2023. Vertical profiles of (a-c) fluorescence spectral backscattering coefficients B λ at 438, 472, 513, 560 and 614 nm and (d-f) aerosol backscattering coefficient β 355 , the ratio B 560 /B 438 and the ratio B 513 /β 355 . Measurements were performed at an angle of 48 • to the horizon.

Figure 10 .

 10 Figure 10. Height profiles of the ratio B 472 /B 438 for smoke episodes during 15 May-20 June 2023 over Moscow. Smoke layers start above 4 km and go up to 11 km.

Figure 11 .

 11 Figure 11. Impact of the aerosol fluorescence on the depolarization ratio in the water vapor Raman channel on the nights of 8-9 and 10-11 June 2023 over Lille. Vertical profiles of (a, c) particle backscattering β 532 , fluorescence backscattering β F and fluorescence capacity G F , as well as (b, d) depolarization ratio δW of the water vapor Raman signal and the water vapor mixing ratio ñW .

Figure 12 .

 12 Figure 12. Relative error in water vapor mixing ratio (WVMR) n W ñW induced by the fluorescence as a function of depolarization ratio δW in the water vapor Raman channel for three values of fluorescence depolarization ratio: δ F = 50 %, 60 %, 70 %. The depolarization ratio of water vapor Raman backscatter in the absence of fluorescence is assumed to be δ W = 2 %.

Figure 13 .

 13 Figure 13. Fluorescence measurements over Lille on the night of 26-27 May and 5-6 June 2023. (a, c) Vertical profiles of the fluorescence backscattering β F , the water vapor mixing ratio ñW , (b, d) the depolarization ratio of the water vapor Raman signal δW , the fluorescence depolarization ratio δ F and parameter η, describing the contribution of the fluorescence to the water vapor channel.

Figure 14 .

 14 Figure 14. Impact of smoke fluorescence on the water vapor measurements. Vertical profiles of the fluorescence backscattering β F , water vapor mixing ratio ñW and bias in water vapor channel n W provided by the fluorescence of smoke for episodes on the nights of (a) 26-27 May, (b) 28-29 May and (c) 16-17 June 2023 for the time interval 21:00-02:30 UTC over Lille. (d) Error n W ñW introduced by smoke fluorescence for the three episodes.

Atmos. Meas. Tech., 17, 1023-1036, 2024 https://doi.org/10.5194/amt-17-1023-2024

https://doi.org/10.5194/amt-17-1023-2024 Atmos. Meas. Tech., 17, 1023-1036, 2024

Acknowledgements. We acknowledge funding from the CaPPA project funded by the Agence Nationale de la Recherche (ANR) through the PIA under contract no. ANR-11-LABX-0005-01, the "Hauts de France" Regional Council (project ECRIN) and the European Regional Development Fund (FEDER). The ESA/QA4EO program is greatly acknowledged for supporting the observation activity at LOA. The work from Qiaoyun Hu was supported by the ANR (ANR-21-ESRE-0013) through the OBS4CLIM project. Development of the fluorescence lidar in Moscow was supported by the Russian Science Foundation (project no. 21-17-00114). The work of Sergey Khaykin was partly supported by the ANR (21-CE01-335 0007-01 PyroStrat project).

Financial support. This research has been supported by the Agence Nationale de la Recherche (grant nos. ANR-11-LABX-0005-01, ANR-21-ESRE-0013 and ANR-21-CE01-0028-02) and the Russian Science Foundation (grant no. 21-17-00114). Publisher's note: the article processing charges for this publication were not paid by a Russian or Belarusian institution.

Review statement. This paper was edited by Simone Lolli and reviewed by Alexandros D. Papayannis and one anonymous referee.

Data availability. Lidar measurements are available upon request (philippe.goloub@univ-lille.fr).

Author contributions. IV processed the data and wrote the paper. QH and TP performed the measurements in Lille. PG supervised the project and helped with paper preparation. WB modified LILAS for polarization measurements. MK and NK performed the measurements in Moscow. SK analyzed transport of smoke layers, and RM derived RH profiles from lidar measurements.

Competing interests. The contact author has declared that none of the authors has any competing interests.

Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


Sunrise Wang
email: sunrise.wang@centralesupelec.frnicolas.gac

Nicolas Gac

Hugo Miomandre
email: hugo.miomandre@insa-rennes.fr

Jean-Francois Nezan
email: jean-francois.nezan@insa-rennes.fr

Karol Desnos
email: karol.desnos@insa-rennes.fr

Francois Orieux
email: francois.orieux@universite-paris-saclay.fr

An Initial Framework for Prototyping Radio-Inteferometric Imaging Pipelines

Keywords: Algorithm Design Space Exploration, SKAO, PREESM, Radio-Interferometry, Resource Estimation

Although large radio-telescope arrays allow us to observe the celestial sphere with an unprecedented level of detail and sensitivity, additional antennas drastically increases the cost of processing and storing their data, complicating the design of computing hardware. Our overall goal is to provide a system, which we term SimSDP, to aid in their design. It will achieve this by providing resource usage estimations for some given imaging pipeline and hardware architecture, allowing for more informed decisions when building the production systems. We lay the groundworks in this paper by presenting and validating an initial system that implements three different imaging pipelines. We find that in most cases, our system is able to accurately estimate the scaling across both algorithmic parameters as well as parallelization when compared to measured data, with errors roughly in the 1-5% range, demonstrating its ability to inform design decisions.

Introduction

Radio-telescopes capture information of our skies within the radio spectrum, allowing for the study of a host of otherwise invisible natural phenomena, with some examples including unionized gas clouds, and low-frequency synchotron radiation. Although many initial radio-telescopes were single-dish instruments, modern advancements in the field of aperture synthesis allows us to instead use antenna arrays with radio-interferometry. This is advantageous as it allows for much larger apertures and higher sensitivities, allowing us to observe smaller and fainter objects.

The currently under construction Square Kilometer Array (SKA) is an example of such a telescope, and on completion, will be the largest antenna array in the world. An issue with the SKA, as well as other large antenna arrays, is their computational cost, as the raw data obtained from the antennas is processed digitally. Thus, the more antennas, the larger the amount of data to be processed, with estimates for the SKA being around 0.4TB/s, leading to roughly 34.5 PBs per day. Storing this data for any significant period of time is prohibitively expensive, leading to strict time restrictions during processing. This, coupled with the fact that the type of science performed dictates the algorithmic pipeline, makes hardware and software design particularly challenging.

To aid in the above, a system that provides resource usage estimations for various radio-interferometric algorithmic pipelines and hardware architectures is desired. Unfortunately, as far as we are aware, none currently exists. The closest to this is DALiuGE [START_REF] Wu | DALiuGE: A graph execution framework for harnessing the astronomical data deluge[END_REF], a graph-based execution framework targeted at radio-astronomy, and DASK [START_REF] Rocklin | Dask: Parallel Computation with Blocked algorithms and Task Scheduling[END_REF], which has been used in radio-interferometry systems such as RASCIL [START_REF] Cornwell | Radio astronomy simulation, calibration and imaging library[END_REF]. However, from our understanding, these are primarily execution frameworks and do not provide easy resource usage estimations, particularly at scale.

It is our aim to introduce such a system, which will be termed SimSDP as it is targeted at the design of the SKA supercomputer, the Science Data Processor (SDP). This system will consist of two main parts. The first is responsible for both the algorithmic descriptions, as well estimating the performance for hardware at a local level, and the second aims to generalize the performance estimations across different hardware architectures.

This paper lays the groundwork for the former, by presenting a generic framework for radio-inteferometric algorithms, and providing three implementations of imaging pipelines within. We evaluate our framework's ability to aid in the design of the SDP by estimating the resource usage of these different pipelines across a variety of parameters and two levels of parallelization, and compare these against measured data to see if we can draw similar conclusions.

The main contributions of this paper are:

-A generic framework describing radio-interferometric algorithms, together with the full data-flow description and implementation of three different imaging pipelines, -An evaluation of our framework's ability to aid in the design of the SDP, by comparing both estimated resource usages against measured across a variety of parameters and two parallelization levels,

The remainder of this paper is structured as follows: Section 2 provides a brief introduction to radio-interferometry as well as the generic pipeline, Section 3 discusses our implementation of this, and provides details on our implemented algorithms and how we perform our estimations, Section 4 discusses our results, and Section 5 concludes our work.

Radio-Interferometric Imaging

Radio-interferometers measure the sky using antenna arrays. Samples are produced by pairs of antennas in the array, termed baselines. Each sample, termed a visibility, is the instrumental response for some given time duration and electromagnetic frequency for a specific baseline. Visibilities can be defined with the Radio-Interferometric Measurement Equation (RIME) [START_REF] Smirnov | Revisiting the radio interferometer measurement equation: I. A full-sky Jones formalism[END_REF]:

V (u, v) = G uv D uv (l, m) 1 √ 1-l 2 -m 2 I(l, m)e -2πi(ul+vm+w( √ 1-l 2 -m 2 -1)) dldm (1)
where G denotes the direction independent effects, such as antenna gain, D denotes the direction dependent effects, such as Faraday rotation caused by the Earth's ionosphere, (u, v) refers to the difference between antenna positions in the frame of the Earth's rotation, (l, m) refers to spatial coordinates on the celestial sphere, and I is the true sky.

If the D and e w( √ 1-l 2 -m 2 -1) terms are ignored, Equation 1 simplifies to

V (u, v) = G uv I(l, m)e -2πi(ul+vm) dldm (2)
a 2-dimensional Fourier transform, allowing us to retrieve an image of the true sky through an inverse Fourier transform. The image produced by the inversion of Equation 2, termed the dirty image, contains artefacts caused both by the partial sampling of the Fourier domain, as antenna arrays are inherently sparse, as well as by the omission of the w and D terms. The former results in a convolution between the Fourier transform of the sampling pattern, ie. the Point Spread Function (PSF), and the real sky, while the latter results in non-isoplanatic image-plane distortions. Radiointerferometric imaging aims primarily to correct these in order to produce an image use-able for scientific purposes.

The Radio-Inteferometric pipeline [START_REF] Clark | An efficient implementation of the algorithm 'CLEAN[END_REF] employs a nested loop structure (Figure 1). The outer ie. major loop ∆ calculates the difference δ în between some estimated sky image în and the measured values v, with î0 typically being a blank image.

Although there are many possible approaches to ∆, the most common general method currently involves transforming în to the same domain as v using a degridding operator, which typically includes a Fourier transform and an extrapolation operator, performing the subtraction to obtain the difference in visibility space, and then using a gridding operator which includes an inverse Fourier transform together with an interpolation operator to obtain δ în .

In addition to computing the difference, the major loop is also responsible for correcting for the w and D terms. This can be achieved either in tandem with the gridding convolution stage G † [START_REF] Cornwell | The Noncoplanar Baselines Effect in Radio Interferometry: The W-Projection Algorithm[END_REF][START_REF] Bhatnagar | Correcting direction-dependent gains in the deconvolution of radio interferometric images[END_REF][START_REF] Van Der Tol | Image Domain Gridding: a fast method for convolutional resampling of visibilities[END_REF] or through various discretization methods [START_REF] Ord | Interferometric Imaging with the 32 Element Murchison Wide-Field Array[END_REF][START_REF] Cornwell | Radio-interferometric imaging of very large fields-the problem of non-coplanar arrays[END_REF][START_REF] Offringa | Wsclean: an implementation of a fast, generic wide-field imager for radio astronomy[END_REF]. There has also been work that employs a hybrid of the two strategies [START_REF] Cornwell | Wide field imaging for the square kilometre array[END_REF][START_REF] Ye | High accuracy wide-field imaging method in radio interferometry[END_REF].

Finally, the major loop also allows for calibration of the direction independent effects G, which we don't show in Figure 1 as it is not the focus of our work.

The inner ie. minor loop Ψ is responsible for removing the convolution artefacts caused by the partial sampling of the Fourier domain from δ în , and then combining this information with în to produce an estimate în+1 for the next major cycle iteration. This is achieved by employing a deconvolution algorithm, with Fig. 1. The general framework for producing the final sky estimate î is iterative, where în is updated while taking into account the difference δ în between it and the measurements v. The degridding-gridding step (∆) is responsible for computing this, whereas the deconvolution step Ψ is responsible for generating the next estimate în+1.

some popular examples being CLEAN [START_REF] Högbom | Aperture synthesis with a non-regular distribution of interferometer baselines[END_REF] and its variants [START_REF] Schwab | Relaxing the isoplanatism assumption in self-calibration; applications to low-frequency radio interferometry[END_REF][START_REF] Cornwell | Multiscale CLEAN Deconvolution of Radio Synthesis Images[END_REF][START_REF] Rau | A multi-scale multi-frequency deconvolution algorithm for synthesis imaging in radio interferometry[END_REF], and convex optimization methods that regularize based on entropy [START_REF] Cornwell | A simple maximum entropy deconvolution algorithm[END_REF] or sparsity [START_REF] Wiaux | Compressed sensing imaging techniques for radio interferometry[END_REF][START_REF] Carrillo | PURIFY: a new approach to radio-interferometric imaging[END_REF][START_REF] Ammanouil | A parallel and automatically tuned algorithm for multispectral image deconvolution[END_REF].

As the deconvolution methods can be partial or contain errors in their estimation, the new estimate în+1 is passed back to the major loop. This is done until δ în is predominantly noise.

Implementation

We implement the generic structure defined in Figure 1, as well as three different imaging algorithms using PREESM [START_REF] Pelcat | Preesm: A dataflow-based rapid prototyping framework for simplifying multicore DSP programming[END_REF]. We find PREESM to be particularly well suited to our needs, as algorithms are defined in data-flow diagrams, which allows for easy definition of additional imaging algorithms. It also allows for automated parallelism, which we use for our different levels of parallelization. Finally, PREESM also allows for easy estimation of resource consumption when provided with the relevant data for each actor.

The generic top-level data-flow diagram provides a common interface for algorithms defined in ∆ and Ψ ie. degridding-gridding and deconvolution. This enables us to encompass a large number of different imaging pipelines, as most work in the domain tend to focus on either.

The rest of this section will look to provide a high-level overview on the algorithms we implemented for both ∆ and Ψ , as well as how parallelization is performed. It also discusses how we perform our resource usage estimations. We provide our system and data in our on-line repository3 .

Algorithms for ∆ and Ψ

We implement three different degridding-gridding algorithms (∆) in our pipeline. The first is the Direct Fourier Transform (DFT) degridder, which transforms the image directly to the visibilities using:

V (u, v) = L,M l=0,m=0 I(l, m)e -2πi(ul+vm) (3)
As it is unnecessary to perform the equation for areas where I is zero, the overall complexity for this method is O(n v n s ), where n v is the total number of visibilities, and n s is the number of sources. Following this is the Fast-Fourier Transform (FFT) degridder, which first performs a Fast Fourier Transform on the image to obtain the gridded visibilities, which are then subsequently convolved onto the continuous visibility positions. This algorithm has a complexity of O(n 2 g log 2 n g + n v n dgk ), where n g is the support size of the grid, and n dgk is the support size of the de-gridding kernel.

Both the DFT and FFT degridders also perform a gridding operation after the subtraction of visibilities, which has a complexity of O(n 2 g log 2 n g + n v n gk ) where n gk is the size of the gridding kernel.

The last implemented degridding-gridding algorithm is the more recent Grid to Grid (G2G) method by Monnier et al. [START_REF] Monnier | Fast Sky to Sky Interpolation for Radio Interferometric Imaging[END_REF]. This method performs degridding and gridding in a single step, and performs the subtraction in image space rather than visibility. The total complexity of this algorithm is O(2n 2 g log 2 n g + 2n v (n gk + n dgk )), the same as the FFT degridder. However, this algorithm is more performant, as it employs a visibility simplification stage which diminishes the size of n v . Furthermore, the subtraction is done in image space, which may be faster as the number of visibilities often dwarfs the number of pixels in the image.

We choose to focus on degridding-gridding (∆) algorithms for this work, thus only implement a single deconvolution algorithm (Ψ ), that being Högbom CLEAN. [START_REF] Högbom | Aperture synthesis with a non-regular distribution of interferometer baselines[END_REF], which we selected for its simplicity. This algorithm iteratively finds the brightest source in the image, and subtracts the PSF from its position, terminating after a preset number of minor cycles n m . The complexity of this algorithm is O(n 2 g +n m n p ), where n p is the support of the PSF. As this algorithm finds a single point source per iteration, n s = n m , which is pertinent for the DFT degridder.

Parallelization

To parallelize ∆, we design the data-flow diagram to divide the visibilities based on some specified number of threads, grid or degrid each subset of visibilities independently, and then merge each thread's results. This manner of parellelization is natural for ∆ as it requires looping through the visibilities. We don't parallelize Ψ as Högbom CLEAN has cross iteration data dependencies, drastically complicating process.

We create two different levels of parallelization with 1 and 4 cores. We set the data division parameter in our data-flow diagrams to not divide the visibilities in the former case, and divide them into 4 disjoint subsets in the latter.

Estimating Resource Utilization

We perform computation time and memory estimations for our system. For the former, PREESM requires time information for each actor, either in the form of an exact value, or some function dependent on the input data. We opt for the latter as it affords us the ability to predict how the algorithms scale.

We perform benchmarking on optimized versions of the code for each actor, while varying the number of visibilities, the image size, and the number of minor cycle iterations. We then fit polynomials of the appropriate degree to these, which we then transfer to PREESM. We average 30 samples for each benchmark data point.

We perform the benchmarking for actors that contain simple and nearoptimal implementations directly in PREESM. For more complex actors, such as the G2G and FFT degridding-gridding algorithms, we profile the code of Monnier et al. [START_REF] Monnier | Fast Sky to Sky Interpolation for Radio Interferometric Imaging[END_REF].

The dataset used to generate these benchmarks is obtained from the SKAO simulated database repository. We provide full information on this in Section 4.1. In order to vary the number of visibilities, we either truncate or duplicate them. We opt for this over random generation because the antenna layout may have an effect on the performance of some algorithms, such as G2G.

The machine used to obtain the benchmarks is a Dell Precision 3551, with an Intel Core i7 10875H 10th generation processor and 64GB of RAM. Re-profiling and fitting is required for other hardware as our fitted models become invalid.

Unlike the computational time estimations, memory estimations in PREESM are much less involved, with PREESM simply reporting the maximum allocated inter-actor memory. It does not estimate intra-actor memory as it does not perform any code analysis, and primarily serves as a lower-bound estimation.

Results and Discussion

This section provides details on the simulated dataset that we use for our experiments, our framework validation, as well as our resource estimation results.

Simulated Database

The simulated dataset used for our validation and experiments is the same as the one used for our benchmarking, and is the small database from the SKA simulated dataset repository4 . It is simulated to have a field of view of 1 degree, with 33 point sources within. The array has a total of 512 antennas, leading to 130816 baselines. The visibilities were sampled from each baseline once every 30 seconds, for a total of 15 minutes, creating 3924480 total visibilities. This dataset is single polarization single wavelength, which is well suited to the current limitations of our system, as it is currently unable to handle different polarizations, extended sources, and multiple frequency channels.

Pipeline Validation

RASCIL DFT

Fig. 2. Final deconvolved images of RASCIL and our DFT degridder-gridder. We can see that the sources appear in the same positions with similar energy ratios, showing that our pipeline functions as intended. Our other images are similar and available on our repository.

We validate our pipeline by performing a visual inspection of its reconstructed images to the ones produced from the well known imaging framework RAS-CIL [START_REF] Cornwell | Radio astronomy simulation, calibration and imaging library[END_REF]. Figure 2 shows the images produced both by RASCIL, and our DFT degridder. The other algorithms produce very similar images, and are available on our repository.

One can see that the general structure of the produced image is the same, with the point sources and their energy ratios being similar. The differences are primarily due to minor algorithmic discrepancies, such as parameters in the deconvolution algorithm.

Computing Time Estimations

We evaluate our computation time estimations by comparing them to the measured times of our pipelines obtained using the GNU time tool. We plot these against different numbers of grid cells, numbers of visibilities, and the numbers of minor cycles per major iteration. We found our estimations to accurately predict the increase in time with increasing algorithm parameters, with most cases having an average error within the 1-5% range, allowing us to draw similar scaling conclusions to the measured data. For example, we can see that the DFT algorithm scales much poorer compared to the others with regards to the number of sources in the image, as well as the number of visibilities.

We also found that our framework allows us to draw valid conclusions with regards to scaling in terms of parallelization, with it accurately predicting in most cases the roughly 1.5 -3× speedup obtained when parallelizing across 4 cores.

Despite the positives, we found various areas in which the accuracy of our estimations was lacking. This can be seen for the FFT and G2G algorithms when varying the number of grid cells, and DFT in the cases where the minor cycles are varied. There are three main reasons for these poor estimates:

-With the exception of the DFT, the code benchmarked was much more optimized than the code measured; -Our fitting functions introduce additional error; -PREESM does not perform static analysis.

Although the former two points can be easily resolved by taking more accurate and a larger number of benchmark samples, the latter is a drawback to our current framework, which only performs static analysis. This is particularly evident in the results for G2G, which performs a visibility simplification step to dynamically reduce the number of visibilities needing to be processed.

Memory Usage Estimations

We evaluate the memory usage estimations by comparing them against the measured maximum resident set size obtained by the GNU time tool in verbose mode. We plot the memory usage against both the number of grid cells and the number of visibilities for both levels of parallelization. Unlike the computational time experiments, we do not vary the number of minor cycles as it does not have a major effect on memory usage.

Figure 4 shows the estimated and measured memory usage in our system. Much like the computation time estimations, our framework predicts well the increase when scaling the amount of data, with an error of 1-5% for all test cases. It also predicts the 1.1 -1.5× increase in memory usage when increasing parallelization.

These results are not surprising, as PREESM is responsible for allocating the inter-actor memory. However, it is unexpected that the measured memory usage is lower than the estimated. As mentioned in Section 3.3, the estimate should provide a lower-bound, which is in conflict with what we observed. We are still investigating possible explanations for this discrepancy, and will leave a more detailed analysis for future work.

Conclusions and Future Work

In this paper, we lay the groundworks for SimSDP, a system aimed at aiding in the design of supercomputers processing data from large radio-telescope arrays by providing resource usage estimations at various scales and hardware architectures. We achieve this by performing a study on a preliminary system aimed at such, which includes three different algorithmic pipelines, and two levels of parallelization. We found that the estimations obtained from our framework allows us to draw similar scaling conclusions to the measured results. This gives impetus in continuing the design and implementation of SimSDP.

Drawbacks

One major drawback of our system is that it currently only performs static analysis, which may lead to it providing poor estimates in cases where algorithms perform dynamic reductions to the data being processed. Also in line with this, there are additional accuracy improvements that can be made to our estimations through more and better benchmarking data. Finally, our system currently does not handle many types of datasets, as it doesn't currently support objects containing extended sources (eg. nebulae), different polarizations, and multiple wavelengths.

Future Work

Much work still needs to be done for the realization of SimSDP. This includes: extending the current system to support different hardware architectures; implementing more algorithms for ∆ and Ψ , which will allow for our system to support more types of datasets; introduce dynamic analysis by employing SPIDER [START_REF] Miomandre | Demonstrating the SPIDER Runtime for Reconfigurable Dataflow Graphs Execution onto a DMA-based Manycore Processor[END_REF]; estimating energy consumption; and comparing our system's estimations to measurements obtained from production systems [START_REF] Cornwell | Radio astronomy simulation, calibration and imaging library[END_REF][START_REF] Offringa | Wsclean: an implementation of a fast, generic wide-field imager for radio astronomy[END_REF][START_REF] Casa | CASA, the Common Astronomy Software Applications for Radio Astronomy[END_REF].

Fig. 3 .

 3 Fig.3. Computational times of both our estimations (dashed lines) and measurements (solid lines) for our implemented methods across both levels of parallelization.

Figure 3

 3 Figure3shows our results, where the dashed line denotes estimated times, the solid line denotes measured, the lighter colours denote the single core setup, and the darker denote the quad core.

Fig. 4 .

 4 Fig. 4. Memory estimations (dashed lines) and measurements (solid lines) for the various algorithms implemented in PREESM plotted against both number of visibilities, and number of grid cells.

https://gitlab.com/igorawratu/an-initial-framework-for-prototyping-radiointeferometric-imaging-pipelines

https://gitlab.com/ska-telescope/sim/sim-datasets

Acknowledgements This work was supported by DARK-ERA (ANR-20-CE46-0001-01)