Aerosol-related effects on the occurrence of heterogeneous ice formation over Lauder, New Zealand ∕ Aotearoa

The presented study investigates the efficiency of heterogeneous ice formation in natural clouds over Lauder, New Zealand / Aotearoa. Aerosol conditions in the middle troposphere above Lauder are subject to huge contrasts. Clean, pristine air masses from Antarctica and the Southern Ocean arrive under southerly flow conditions, while high aerosol loads can occur when air masses are advected from nearby Australia. This study assesses how these contrasts in aerosol load affect the ice formation efficiency in stratiform midlevel clouds in the heterogeneous freezing range (−40 to 0 ^∘C). For this purpose, an 11-year dataset was analyzed from a dual-wavelength polarization lidar system operated by National Institute of Water and Atmospheric Research (NIWA), Taihoro Nukurangi, at Lauder in collaboration with the National Institute for Environmental Studies in Japan and the Meteorological Research Institute of the Japan Meteorological Agency. These data were used to investigate the efficiency of heterogeneous ice formation in clouds over the site as a function of cloud-top temperature as in previous studies at other locations. The Lauder cloud dataset was put into context with lidar studies from contrasting regions such as Germany and southern Chile. The ice formation efficiency found at Lauder is lower than in polluted midlatitudes (i.e., Germany) but higher than, for example, in southern Chile. Both Lauder and southern Chile are subject to generally low free-tropospheric aerosol loads, which suggests that the low ice formation efficiency at these two sites is related to low ice-nucleating-particle (INP) concentrations. However, Lauder sees episodes of continental aerosol, more than southern Chile does, which seems to lead to the moderately increased ice formation efficiency. Trajectory-based tools and aerosol model reanalyses are used to relate this cloud dataset to the aerosol load and the air mass sources. Both analyses point clearly to higher ice formation efficiency for clouds which are more strongly influenced by continental aerosol and to lower ice formation efficiency for clouds which are more influenced by Antarctic/marine aerosol and air masses.

1 Introduction

Comprehension of cloud formation and cloud lifecycle is essential for adequately describing the Earth's climate system (Stephens, 2005). An important step towards such understanding lies in disentangling and eventually quantifying the complex processes of aerosol–cloud interaction (Lohmann and Feichter, 2005). Aerosol particles are required in the formation of liquid droplets and ice crystals from vapor phase by acting as cloud condensation nuclei (CCN) and as ice-nucleating particles (INPs) in the heterogeneous freezing temperature range from approximately −40 to 0 ^∘C, respectively. In this temperature range, mixed-phase clouds can exist. They are of special interest in current atmospheric research because they can lead to both cooling or warming at Earth's surface (Sun and Shine, 1994; Hogan et al., 2003; Shupe and Intrieri, 2004; Zuidema et al., 2005) and because of their critical role for precipitation (Korolev et al., 2017).

In observational studies, it is difficult to separate thermodynamical, dynamical, and aerosol-related effects on cloud microphysics. Dynamics and thermodynamics are considered to preside over aerosol-related effects on cloud microphysical properties. Nevertheless, significant influence of aerosol conditions on cloud properties, cloud evolution, and precipitation has been confirmed in observational and modeling studies (Seifert et al., 2012; Possner et al., 2017; Solomon et al., 2018; Zhang et al., 2018; Desai et al., 2019).

The strong contrast in aerosol particle load and composition between the Southern Hemisphere and Northern Hemisphere midlatitudes (Tegen et al., 1997) can therefore be expected to lead, as well, to differences in cloud abundance and cloud properties in those regions. Indeed, the recurrent abundance of supercooled liquid cloud layers over the Southern Ocean (Hu et al., 2010; Tan et al., 2014; Huang et al., 2015) has a sensitive impact on the energy balance of this area (Trenberth and Fasullo, 2010; Bodas-Salcedo et al., 2014, 2016, 2019). These regional differences in the occurrence of supercooled liquid clouds might be caused by the lack of INPs in the largely pristine Southern Ocean (Hamilton et al., 2014). Shipborne INP measurements presented by Welti et al. (2020) and McCluskey et al. (2018a) confirm the rareness of INPs in the Southern Ocean. Tatzelt et al. (2022), also from shipborne INP measurements, conclude based on correlation between aerosol physicochemical properties and proxies of biological activity in the ocean as well as backward-trajectory analyses that indications for local INP sources are rare. Therefore, Tatzelt et al. (2022) interpret the INP populations as mixed long-range-transported biogenic and mineral dust INP from coastal or terrestrial sources which are sparse and distant (Vergara-Temprado et al., 2017).

On the other hand, marine organic aerosol is considered to be potentially an important source of INPs, especially over remote oceanic regions such as the Southern Ocean (Wilson et al., 2015; DeMott et al., 2016). By introducing marine organic aerosol into an Earth system model via a physically based parameterization instead of empirical proxies of biological activity, Zhao et al. (2021) concluded that marine INP dominate primary ice nucleation in the Southern Ocean at pressure levels above 400 hPa. McCluskey et al. (2019) compared shipborne INP measurements with model representations and found sea spray aerosol as the dominant contributor to primary ice nucleation in the Southern Ocean, though, also stressing the potential role of long-range-transported mineral dust and the need for modeling and observational studies on vertical and seasonal variability in INP and aerosol composition over remote regions.

The representation of the abovementioned contrasts in weather and climate models may impact these model's performance. Coupled Model Intercomparison Project (CMIP) models consistently showed and still show a positive bias in the absorbed shortwave radiation over the Southern Ocean (Trenberth and Fasullo, 2010; Grise et al., 2015; Zelinka et al., 2020; Varma et al., 2020; Kajtar et al., 2021; Schuddeboom and McDonald, 2021). This is largely attributed to simulating too much ice in Southern Ocean clouds compared to observations, which is another common bias in CMIP models (Kuma et al., 2023). A reduced frequency of liquid clouds leads to too less reflected shortwave radiation. This leads to an enhanced heating at the surface, which in turn causes a too large ocean heat uptake (Franklin et al., 2013; Hyder et al., 2018). The reasons for these observed differences and widespread model biases remain partly unresolved and therefore disputed, but hemispheric contrasts in aerosol load and type, including the low abundance of INPs from terrestrial sources, in combination with atmospheric dynamics, are suspected to play a role (Radenz et al., 2021a).

Mallet et al. (2023) and Humphries et al. (2023) highlight the heterogeneity of the pristine Southern Ocean region and suggest coordinated, long-term, and contrasting observational studies for validating and improving the models. Already, Minikin et al. (2003) performed aircraft in situ measurements over the Northern Atlantic and Southern Ocean and found a factor of 2–3 lower aerosol concentration in the free troposphere over the latter compared to over the former. Gayet et al. (2004) used the same measurements, focusing on cirrus evolution and properties, and derived lower cirrus ice crystal number concentrations in the southern than in the northern midlatitudes. Using spaceborne polarization lidar, Choi et al. (2010) demonstrated a negative correlation of supercooled cloud fraction with dust aerosol frequency. However, Mace et al. (2021a) used ship-based lidar observations (Mace et al., 2021b; McFarquhar et al., 2021) to underline that such spaceborne lidar phase statistics of mixed-phase clouds have a low bias in mixed-phase cloud occurrence over the Southern Ocean. To try to overcome spaceborne lidar observation and detection limits, Villanueva et al. (2020) used spaceborne observations in combination with aerosol reanalysis and found cloud glaciation to be differently susceptible to otherwise similar dust loads in the Southern Hemisphere compared to the Northern Hemisphere, which they found to agree with other findings on hemispheric differences in the mineralogical composition and related freezing efficiency of dust.

In the Southern Hemisphere midlatitudes, ground-based lidar was first used by Kanitz et al. (2011) to assess the thermodynamic phase of stratiform mixed-phase clouds observed over Punta Arenas, Chile. They found 50 % of clouds glaciated at −30 ^∘C, compared to the same glaciation occurring already at −13 ^∘C in the Northern Hemisphere midlatitudes (Seifert et al., 2010). Alexander and Protat (2018) corroborated these results in a similar study performed with lidar data from Kennaook / Cape Grim. Radenz et al. (2021a) performed a contrasting study using a longer-term dataset (than Kanitz et al., 2011) again from clean and pristine Punta Arenas but also including Doppler wind lidar and Doppler cloud radar, from dust-laden Eastern Mediterranean (Limassol, Cyprus), and from polluted northern midlatitudes (Leipzig, Germany). They found the ice formation efficiency for clouds observed over Punta Arenas to be most strongly suppressed compared to Leipzig and Limassol for clouds decoupled from the surface aerosol reservoir. Therefore, Radenz et al. (2021a) emphasize, besides the impact of gravity waves, the importance of quantifying terrestrial INP sources in the Southern Ocean. These terrestrial emissions need to be quantified by observations of free tropospheric aerosol up- and downwind of major landmasses in order to identify regions where cloud microphysics are most strongly impacted by these INPs (Radenz et al., 2021a).

This is strong motivation to extend such investigations of aerosol effects on mixed-phase clouds to New Zealand / Aotearoa. There, besides clean marine conditions, episodic aerosol transport from Australia (wildfire smoke, desert dust, urban haze) is occurring (Liley et al., 2001; Jones et al., 2001; Liley and Forgan, 2009). This enables a unique contrasting study at a single location which would not be possible in the permanently polluted Northern Hemisphere midlatitudes or in the permanently clean environment (of the free troposphere) of, for example, Punta Arenas. To investigate the efficiency of heterogeneous ice formation in natural clouds, this study makes use of a dataset which is unparalleled in the Southern Hemisphere, consisting of long-term ground-based polarization lidar measurements on the South Island /Te Waipounamu of New Zealand.

The paper is structured as follows. In Sect. 2, the used data and methods are introduced. In Sect. 3, first, a measurement example is described as a case study (Sect. 3.1). Later, general air mass statistics based on trajectory and aerosol model reanalysis are provided (Sect. 3.2). The retrieved cloud statistics are then presented in Sect. 3.3. In Sect. 3.4, the separation of clouds based on the air mass statistics is presented. Finally, in Sect. 4, the findings are summarized and an outlook is provided.

2 Data and methods

2.1 Lauder lidar system

Since 1992, lidar observations have been conducted by the National Institute of Water and Atmospheric Research (NIWA), Taihoro Nukurangi, at the Lauder Atmospheric Research Station on the South Island / Te Waipounamu of New Zealand, in collaboration with the Meteorological Research Institute of the Japan Meteorological Agency (MRI/JMA) (Uchino et al., 1995). In February 2009, the lidar system was upgraded with polarization capability (Nagai et al., 2010). After this upgrade, lidar observations are conducted by NIWA in collaboration with the National Institute for Environmental Studies (NIES) in Japan and MRI/JMA. The detailed specifications of the zenith-pointing dual-wavelength (1064, 532 nm) polarization lidar system are described in Nakamae et al. (2014) and Sakai et al. (2016). It uses a Nd:YAG (neodymium-doped yttrium aluminum garnet) laser with a pulse energy of 150 mJ, a pulse width of 5 ns, and a repetition rate of 10 Hz. The telescope diameter is 30.5 cm, and the range resolution is 7.5 m. The field of view is 1 mrad, and the beam divergence is 0.2 mrad. The signal at 1064 nm wavelength is detected with an analogue-mode avalanche photodiode, and two polarization components of the signals at 532 nm wavelength are detected with three photomultiplier tubes. This allows for the retrieval of the particle backscatter coefficient at 532 nm and 1064 nm wavelength using the Klett–Fernald–Sasano algorithm (Klett, 1981, 1985; Fernald, 1984; Sasano et al., 1985), as well as the particle depolarization ratio at 532 nm wavelength. Lauder Atmospheric Research Station is located at 370 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$; 45.037^∘ S, 169.683^∘ E. For this study, an 11-year polarization lidar dataset from 15 March 2009 to 5 April 2020, particularly the co- and cross-polarized signals at 532 nm wavelength, was used.

2.2 Identification of cloud layers and their phase

The cloud layers observed by polarization lidar were classified as either ice-containing or pure liquid. The used methodology was previously applied (Ansmann et al., 2009; Seifert et al., 2010, 2011, 2015; Kanitz et al., 2011) and is therefore only briefly outlined here. A cloud is defined as individual when it is separated from another cloud by more than 5 min in time and 500 m vertically in space. If the distance in space and time is smaller, the cloud layers are counted as only one. The cloud boundaries (start and end times, cloud base, and top height) are determined by visual inspection of the range-corrected lidar signal with the aid of the volume depolarization ratio. Since the lidar signal might be already attenuated at the apparent cloud top, a quality flag (“well-defined cloud top”) is introduced based on the signal strength above the visually defined cloud top where the signal-to-noise ratio must at least be 3 % larger than the background noise level. Next, the volume depolarization ratio profile of each cloud layer is examined to identify ice-containing cloud layers. Irregularly shaped ice crystals depolarize light due to internal reflections. Liquid cloud droplets are spherical and do not depolarize light as long as no multiple scattering (at multiple droplets) occurs (Jimenez et al., 2020). A highly accurate determination of the volume linear depolarization ratio is not necessary for this study. A clear discrimination of liquid water droplet depolarization (typically < 0.2, including multiple scattering effects) and ice crystal depolarization (usually > 0.4) is sufficient (Ansmann et al., 2009). Specular reflection by aligned and falling ice crystals can occur due to the zenith-pointing laser beam (Platt, 1978). This can lead to volume depolarization ratios close to zero, which are not anymore distinguishable from depolarization ratios resulting from backscattering by spherical liquid droplets. These aggravating effects (multiple scattering, specular reflection) need to be kept in mind during cloud phase determination with polarization lidar. Therefore, the depolarization ratio values at the cloud base need to be considered when classifying cloud phase based on depolarization ratio (Seifert et al., 2010; Westbrook et al., 2010; Uchino et al., 1988). Other clear indications of ice-containing clouds are ice virgae which can be unambiguously identified by ground-based lidar (Kanitz et al., 2011).

The efficiency of heterogeneous ice nucleation is mainly controlled by temperature (DeMott et al., 2015). Therefore, the coldest part of a cloud (i.e., the cloud top) offers the most favorable conditions for the start of heterogeneous ice formation (Rauber and Tokay, 1991; Lebo et al., 2008). We thus associate the cloud-top temperature to each of the observed cloud events. In a last step, each categorized cloud gets a cloud-top temperature assigned, obtained from profiles from the GDAS (Global Data Assimilation System) with 1^∘ spatial and 3 h temporal resolution from the National Weather Service's National Centers for Environmental Prediction (NCEP) at the coordinates 45^∘ S, 169^∘ E (GDAS, 2023).

2.3 Auxiliary tools

The Copernicus Atmospheric Monitoring Service – Monitoring Atmospheric Composition and Climate (CAMS-MACC; Morcrette et al., 2009; Inness et al., 2019) model aerosol reanalyses were used to obtain simulated aerosol conditions for each classified cloud. The 3-hourly resolved mass mixing ratios of sea salt aerosol in three particle size ranges (0.03–0.5, 0.5–5, and 5–20 µm radius), mineral dust aerosol in three particle size ranges (0.03–0.55, 0.55–0.9, and 0.9–20 µm radius), hydrophilic and hydrophobic organic matter aerosol, hydrophilic and hydrophobic black carbon aerosol, and sulfate aerosol, as well as air temperature on a 3 × 3 × 12 grid, were used (three longitudes (44.5, 45.25, 46^∘ S), three latitudes (168.5, 169.25, 170^∘ E), and 12 pressure levels (250–1000 hPa), generated using CAMS, 2023). From that grid, each cloud was assigned the corresponding aerosol information (mixing ratios) by trilinear interpolation to the cloud start time and cloud top height above Lauder. Interpolation was carried out in the pressure space. Later, heights were calculated using the barometric formula based on the US Standard Atmosphere (standard temperature of 288.15 K, standard pressure of 1013.25 hPa, and a temperature gradient of 0.65 K per 100 m; U.S. Government Printing Office, 1976).

The HYSPLIT model (HYbrid Single-Particle Lagrangian Integrated Trajectory; Stein et al., 2015; Rolph et al., 2017; HYSPLIT, 2023) was run to calculate backward trajectories for the presented example case and the clustering analysis which is shown in Sect. 3.4. The HYSPLIT backward trajectories were calculated with a starting time at each cloud's starting time rounded to a full hour. The arrival heights above the measurement site were set to the cloud middle heights. After clustering of the backward trajectories using the standard HYSPLIT (2023) clustering tool, each cloud was assigned to its corresponding cluster.

The temporally and vertically resolved air mass source attribution tool TRACE (TRACE, 2023; Radenz et al., 2021b) allows for the retrieval of normed residence times of trajectories within defined height ranges (e.g., the planetary boundary layer or its proxies; the so-called reception height) over certain surface-cover classifications, named geographical areas, and latitude bands. The residence times at each time and height step are summed for each of these classes, where the air parcel was below the reception height. The residence time is an indication for the aerosol characteristics of the particular air parcel. TRACE was run with a temporal resolution of 3 h and the standard reception height of 2 $\mathrm{km}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ for the period March 2009–April 2020 based on HYSPLIT backward trajectories. The following surface-cover classifications (“Water”, “Forest”, “Savanna/shrubland”, “Grass/cropland”, “Urban”, ”Snow/ice”, and “Barren”), latitude bands (“< −60”, “−60 to −30”, “−30 to 0”, “0 to 30”, “30 to 60”, and “> 60^∘”), and named geographical areas (“South America”, “Africa”, “Australia”, “New Zealand”, and “Antarctica/Southern Ocean” (equal to the latitude band “< −60^∘”)) are defined in the TRACE tool (Radenz et al., 2021b). The TRACE normed residence times were linked with the cloud dataset using a nearest-neighbor approach (i.e., closest in time and height to the cloud start times and cloud top heights).

Figure 1 Measurement example from Lauder on 19–20 January 2020 corresponding to the Australian cluster 1 from Fig. 9. Temporal development of the uncalibrated attenuated backscatter (a) and particle depolarization ratio (b) at 532 nm wavelength from 19 January 2020, 12:01 UTC, to 20 January 2020, 11:59 UTC, with indicated liquid (green rectangles) and ice-containing clouds (magenta rectangles). Vertical profiles of the particle backscatter coefficient at 532 nm wavelength (20 January 2020, 04:30–04:40 UTC, light green, evaluated with the Klett–Fernald–Sasano algorithm) (ci); INP concentration obtained with conversion parameters from Mamouri and Ansmann (2016) and the parameterizations of McCluskey et al. (2018b, MC18, marine, light blue) and DeMott et al. (2010, DM10, continental, dark green) (cii). The dashed lines indicate the temperatures from GDAS (20 January 2020, 03:00 UTC). Error bars and shaded areas indicate the uncertainty in the retrieved values. HYSPLIT backward trajectories for 120 h arriving at 3, 5, and 7 km height above Lauder at 19 January 2020, 15:00 UTC (d).

3 Results

3.1 Case study

An example measurement of clouds above Lauder on 19–20 January 2020 UTC is shown in Fig. 1. The temporal developments of the uncalibrated attenuated backscatter coefficient (Fig. 1a) and the particle depolarization ratio (Fig. 1b) at 532 nm show a scene of ice-containing and liquid clouds between roughly 5–9 km height. Up to 3 km height, slightly depolarizing aerosol and at 3.5 km height a lofted aerosol layer were present. Above 9 km height, cirrus clouds were observed. This measurement falls into the period when, above Lauder, elevated smoke layers from the extreme Australian bush fires 2019/20 (Ohneiser et al., 2020; Bègue et al., 2023) were frequently detected. Below 2 km height, some liquid clouds with cloud-top temperatures above 0 ^∘C were measured. Furthermore, low clouds/fog impaired the lidar measurements for a short period on 19 January 2020. During this scene, nine ice-containing clouds (indicated with magenta rectangles) and seven liquid clouds (green rectangles, Fig. 1a and b) were classified.

Figure 1c shows vertical profiles of particle backscatter coefficient at 532 nm wavelength (Fig. 1ci) and INP concentration of marine and continental aerosol (Fig. 1cii). The profile was evaluated for a period without lower-level clouds (20 January 2020, 04:30–04:40 UTC). The INP concentration was estimated from the backscatter coefficient using the method and conversion parameters described in Mamouri and Ansmann (2016) and the freezing parameterizations for continental (DeMott et al., 2010) and marine aerosol (McCluskey et al., 2018b), which were applied for the temperature range from −40 to 0 ^∘C as in Mamouri and Ansmann (2016). For the liquid cloud on 20 January 2020, 04:43–05:01 UTC, with a cloud top at 7.7 km and a cloud-top temperature at −24 ^∘C, this calculation yields INP concentrations at the cloud top of 0.3 L⁻¹ (continental aerosol) and 2 × 10⁻³ L⁻¹ (marine aerosol). For the ice-containing cloud on 20 January 2020, 05:14–12:01 UTC, with a cloud top at 8.9 km and a cloud-top temperature of −34 ^∘C, where freezing is occurring probably regardless of the specific aerosol population, this yields INP concentrations at cloud top of 1.5 L⁻¹ (continental) and 0.2 L⁻¹ (marine). In general, the INP concentration estimates for continental (and, respectively, marine) aerosol have the following range: 2 × 10⁻³ L⁻¹ (4 × 10⁻⁷ L⁻¹) at −10 ^∘C, 0.06 L⁻¹ (1 × 10⁻⁴ L⁻¹) at −20 ^∘C, 0.9 L⁻¹ (0.04 L⁻¹) at −30 ^∘C, and 5 L⁻¹ (5 L⁻¹) at −40 ^∘C. These differences of up to 4 orders of magnitude emphasize the strong impact that the continental aerosol contribution has on INP concentrations and the potential resulting changes in freezing efficiency. For comparison, Radenz et al. (2021a) reported lidar-derived INP concentrations for Punta Arenas of about 3 × 10⁻⁵ L⁻¹ (7 × 10⁻³ L⁻¹) at −20 ^∘C (−30 ^∘C) when only pristine marine aerosol is present, while already small quantities (1 %–4 %) of continental aerosol lead to a significant increase of the INP concentrations to about 0.06 L⁻¹ (0.06 L⁻¹) at −20 ^∘C (−30 ^∘C).

HYSPLIT backward trajectories for 120 h arriving above Lauder at 3, 5, and 7 km height on 19 January 2020, 15:00 UTC, are shown in Fig. 1d. While the lower backward trajectory, probably representative of the near-ground as well as lofted aerosol layer, points towards rather local or regional air mass origin, the upper backward trajectories come from Australia or its direction. Therefore, this measurement example (i.e., the ice-containing cloud on 19 January 2020, 15:59–18:51 UTC, also indicated with one of the magenta rectangles in Fig. 1a and b) belongs to the Australian backward trajectory cluster 1 (Fig. 9a), exhibiting higher ice formation efficiency (Sect. 3.4, Fig. 9b). At lower levels, the air masses were transported slowly due to a weak-gradient high-pressure situation, and only in the upper half of the troposphere (> 6 km height) did a faster flow regime from Australia prevail (NIWA, 2023).

Figure 2 Vertical profiles at 19 January 2020, 15:00 UTC, of TRACE normed residence times of surface-cover (a) and named geographical area classifications (b), MACC aerosol mixing ratios excluding sea salt (c), and temperature and relative humidity from GDAS (d).

Figure 2 shows vertical profiles of TRACE surface-cover and named geographical area classifications, MACC aerosol mixing ratios excluding sea salt, and temperature and relative humidity from GDAS temporally closest (19 January 2020, 15:00 UTC) to the abovementioned ice-containing cloud. The TRACE surface-cover classifications indicate air masses influenced by Water and Grass/cropland in the lower altitudes and by Savanna/shrubland (besides Water) in the higher altitudes (Fig. 2a). The TRACE named geographical area classifications (Fig. 2b) indicate air masses influenced by Australia above 5.5 km height and New Zealand and Antarctica/Southern Ocean up to about 2–2.5 km. Although the backward trajectory arriving at 7 km height of the single HYSPLIT run in Fig. 1d is staying exclusively in the free troposphere, the statistical results of TRACE (having a reception height of 2 $\mathrm{km}\mathrm{a}.\mathrm{g}.\mathrm{l}.$) clearly indicate a continental contribution at these heights (Fig. 2a and b). The MACC aerosol mixing ratios (excluding sea salt) are dominated by organic matter and sulfate which are decreasing above 5.5 km, while dust and black carbon are marginal throughout the column (Fig. 2c). According to the GDAS profile, the temperature at 7 km is about −20 ^∘C at high relative humidity (Fig. 2d).

Figure 3 Seasonal occurrence of the heterogeneous freezing regime at Invercargill / Waihōpai based on radiosoundings. Panel (a) shows the average heights of the 0 ^∘C isotherm (red) and −35 ^∘C isotherm (blue). Frequencies of dew point spreads of 0–2 K (green) and 2–6 K (orange) within the heterogeneous freezing temperature range are shown in panel (b). The number of observations are given on the top of each plot. Soundings from 1950–2022 were used.

3.2 Air mass statistics

Before the presentation of the mixed-phase cloud statistics, general information on the conditions for heterogeneous freezing and the aerosol properties are provided in this section. Long-term records of radiosoundings at Invercargill Airport about 200 km southwest of Lauder were used to assess the abundance of thermodynamic conditions under which heterogeneous freezing occurs, namely, temperatures between 0 and −35 ^∘C and the availability of moisture. Given the short distance to Lauder, we assume that the conditions are comparable to those from Invercargill / Waihōpai, for which the valuable radiosonde dataset is available. Tradowsky et al. (2018) compared radiosonde and automatic weather station data collected at Lauder and Invercargill / Waihōpai and concluded that especially temperature anomalies are well correlated between these two sites. Figure 3 shows the seasonal occurrence of that heterogeneous freezing regime over Invercargill / Waihōpai. The heights of the isotherms show a maximum in the summer season (austral, DJF) with 2.6 km for the 0 ^∘C isotherm and 8.4 km for the −35 ^∘C isotherm (Fig. 3a). In austral autumn (MAM), their heights are only slightly lower, while in winter their heights decrease to 7 km for the −35 ^∘C isotherm and to 1.4 km for the 0 ^∘C isotherm. During spring (SON), the isotherms remain at lower heights than in autumn. Large-scale dew point spreads less than 2 K are usually needed for clouds to form (Poore et al., 1995). Depending on the vertical air motions, dew point spreads up to 6 K could be sufficient as well. Figure 3b shows the relative frequency of dew point spreads between 0–2 K and between 2–6 K, respectively, for heights with temperatures in the heterogeneous freezing regime. The frequencies of the 0–2 K dew point spreads during summer and autumn are nearly equal with 14 %, and the maximum occurs in winter with a frequency of 17 % being only slightly higher than the one in spring. Therefore, suitable conditions for clouds to form in the heterogeneous freezing regime are present throughout the year at Invercargill / Waihōpai. Climatologically, one can assume similar conditions for the site of Lauder, which is just 200 km north of the radiosonde station in Invercargill / Waihōpai. However, due to the vicinity of Lauder to the mountains and associated lee effects, the free troposphere over Lauder might contain less water vapor. Therefore, the average dew point spread might be lower than over Invercargill / Waihōpai. In turn, the likelihood of strong upward motion due to orographic mountain waves might be pronounced over Lauder, which would favor the formation of free-tropospheric cloud layers.

Figure 4 Time series of 28 d mean CAMS-MACC sea salt (blue), dust (orange), organic matter (green), black carbon (black), and sulfate (red) aerosol mixing ratios at 500 hPa above Lauder from January 2009 to April 2020.

Figure 5 Same as Fig. 4 but excluding sea salt.

The general aerosol conditions above Lauder were derived based on the MACC dataset. The MACC mixing ratios of the three size ranges of sea salt and dust aerosol species and the species of organic matter and black carbon were summed up and averaged individually on a 28 d time window over the measurement period from January 2009 to April 2020. Figure 4 shows a time series of these 28 d means of the MACC aerosol mixing ratios at the 500 hPa level above Lauder. Obviously, the simulated aerosol mixing ratios are dominated by sea salt which shows weakly pronounced minima during (austral) summer and maxima during (austral) spring. Figure 5 shows the same data at the same pressure level but excluding the dominant sea salt mixing ratio. In strong contrast to the total aerosol mixing ratios, the remaining aerosol mixing ratios show a distinct annual cycle with maxima during (austral) spring or summer. The most abundant species are organic matter followed by dust and sulfate with marginal amounts of black carbon. The simulated aerosol mixing ratios reach exceptionally high values in austral summer 2019/20 due to the Australian bush fires at that time. Based on long-term sun photometer observations, Liley and Forgan (2009) showed that the median aerosol optical thickness (AOT) above Lauder is about 0.03 at 500 nm wavelength with a minimum during austral winter (0.02) and a maximum during austral spring and summer (0.04). They also report no significant variation of AOT with surface humidity below 90 %. The peak in AOT in spring has been attributed to long-range transport of biomass burning products, as implied by the correlation with carbon monoxide measurements presented by Liley et al. (2001) and Jones et al. (2001). With regard to the usage of the MACC mixing ratios (as an indication of aerosol composition and abundance for the analysis of aerosol impact on the formation of mixed-phase clouds), this dominance of sea salt could be obstructive, and, as the annual cycle in Fig. 4 shows, large amounts of sea salt might rather be a tracer of air containing less other aerosol. Therefore, sea salt is omitted in the subsequent analyses and discussions.

Figure 6 Profiles of normed residence times of trajectories above named geographical areas (a, b) and surface-cover classifications (c, d) in terms of kilometers above ground level ($\mathrm{km}\mathrm{a}.\mathrm{g}.\mathrm{l}.$) (a, c) and temperature (b, d) from the temporally and vertically resolved air mass source attribution tool TRACE (Radenz et al., 2021b).

Results of the temporally and vertically resolved air mass source attribution tool TRACE for trajectories arriving above Lauder are shown in Fig. 6. The height profiles of normed residence times of trajectories above the named geographical areas indicate New Zealand and Antarctica/Southern Ocean as the most important air mass sources up to heights of about 2 km (Fig. 6a). Above that, Australia is the major source of air masses, while Africa and South America show residence times that are an order of magnitude lower. The temperature profiles of named geographical area classifications indicate New Zealand as the major source up to shortly above 0 ^∘C, at lower temperatures Antarctica/Southern Ocean and again Australia become more important (Fig. 6b). Regarding the height profiles of normed residence times of trajectories above surface-cover classifications, Water is dominant over the whole column followed by Savanna/shrubland, except in the lowest 1–2 km where the other categories (except Barren) show slightly higher normed residence times (Fig. 6c). In terms of temperature profiles of surface-cover classifications, a very similar picture as in the height profiles is drawn apart from relatively larger residence times above Snow/ice at about −15 to 0 ^∘C (Fig. 6d). It can be concluded that Lauder is seeing air masses influenced by New Zealand / Aotearoa and Australia, mainly at heights (temperatures) above (below) 2 km (−5 ^∘C), as well as air masses of marine or Southern Ocean origin.

Figure 7 Classified liquid (blue) and ice-containing clouds (red) in the height range 0–15 km a.g.l. above Lauder during the analyzed years 2009–2020 (top to bottom).

Table 1 Overview of the extent of the 11-year Lauder cloud dataset. Number (N) of categorized liquid (liq.) and ice-containing (ice) clouds, including relative time of cloud occurrence (t_sampled) for the full dataset and for the subdataset covering the heterogeneous freezing range −40 to 0 ^∘C and for clouds with a well-defined (WD) cloud top (see Sect. 2).

3.3 Cloud statistics

In total, the 11-year dataset covers more than 9000 h of polarization lidar measurement data. Figure 7 shows all classified cloud layers from 2009 to 2020. As can be seen from the increased number of classified clouds, three periods in 2011, 2015, and from the end of 2019 until 2020 experienced more frequent lidar observations. These were conducted to characterize the plumes of the erupted Chilean volcanoes Puyehue-Cordón Caulle (2011) and Calbuco (2015) and to obtain a larger database specifically for this study (2019/20). A total of 3360 clouds were identified with a total occurrence time of more than 6000 h, with 1576 out of them in the temperature range between −40 and 0 ^∘C. The number of clouds (total, well-defined, in the range of −40 to 0 ^∘C, respectively) as well as the relative time covered by the occurrence time of the identified clouds for each year is listed in Table 1. Overall, about 45 % of the observed clouds (1497) did not have a well-defined cloud top due to too large cloud optical thickness. The remaining 1863 clouds were classified as well-defined (WD). In the cloud-top temperature range of −40 to 0 ^∘C, 61 % of clouds (973) have a well-defined cloud top. The fraction of ice-containing clouds at a specific cloud-top temperature can be considered characteristic for a certain geographical region and is widely used to describe ice nucleation in models or observations (e.g., Seifert et al., 2015; McCoy et al., 2015; Alexander and Protat, 2018; Villanueva et al., 2021; Radenz et al., 2021a). Figure 8 shows the fraction of ice-containing clouds in eight cloud-top temperature intervals of 5 K from −40 to 0 ^∘C for Lauder (this study), Leipzig (Seifert et al., 2010), and Punta Arenas (Kanitz et al., 2011). The error bars denote the statistical uncertainty by the standard error σ defined by the following formula:

$$\begin{array}{}\text{(1)}& \mathit{\sigma }=\sqrt{f(\mathrm{1}-f)/n},\end{array}$$

with f as the frequency of occurrence of ice-containing clouds and n as the total number of observed clouds in the specific temperature range. The fraction of ice-containing clouds including clouds without well-defined cloud tops is slightly larger at most cloud-top temperature ranges and is shown too (Fig. 8, dashed line). Leipzig (51.4^∘ N, 12.4^∘ E) represents typical polluted aerosol conditions of the northern midlatitudes (Seifert et al., 2010; Kanitz et al., 2011), while Punta Arenas (53.1^∘ S, 70.9^∘ W; southern Chile) represents clean and pristine southern midlatitudes (Jimenez et al., 2020; Radenz et al., 2021a). For Lauder, an increase of the relative number of ice-containing clouds is observed at about a temperature range of −20 to −15 ^∘C. In this temperature range, the fraction of ice-containing clouds at Lauder is about a factor of 5.4 smaller than at Leipzig and close to that at Punta Arenas. On the other hand, in the temperature range of −30 to −25 ^∘C, the fraction of ice-containing clouds at Lauder is 1.3 times smaller than at Leipzig and still 1.6 times larger than at Punta Arenas. Lauder and Punta Arenas are subject to generally low free-tropospheric aerosol loads, which suggests that the lower ice formation efficiency at these two sites compared to Leipzig is related to low INP concentrations.

Figure 8 Fraction of ice-containing clouds as a function of cloud-top temperature in intervals of 5 K for Leipzig (red; Seifert et al., 2010), Punta Arenas (blue; Kanitz et al., 2011), and Lauder considering only well-defined clouds (black, this study) and all clouds (grey, this study). The error bars represent the statistical significance calculated with Eq. (1).

Figure 9 Clusters of 1544 (72 h HYSPLIT) backward trajectories arriving at cloud middle height and cloud start time rounded to full hour (the four clusters explain about 85 % of the clustered backward trajectories' total spatial variance) (a); same as Fig. 8 but all well-defined clouds (black, same as in Fig. 8) and well-defined clouds corresponding to cluster 1 (red), cluster 2 (blue), and clusters 3 and 4 (green) (b).

3.4 Separation of clouds based on air mass statistics

Backward trajectory tools were used to investigate the aerosol effect on ice formation efficiency in clouds above Lauder. In a first approach, backward trajectories for each classified well-defined cloud were calculated and clustered (see Sect. 2.3). The resulting clusters are shown in Fig. 9a. Cluster 1 (35 % of trajectories) can be described as direct Australian transport. Contrasting to this, cluster 2 (26 % of trajectories) represents trajectories with fast transport of Antarctic air masses. Clusters 3 (18 %) and 4 (21 % of trajectories) contain trajectories of air masses transported above ocean and, therefore, can be characterized as marine. Figure 9b shows the fraction of ice-containing clouds against cloud-top temperature for all well-defined clouds and for well-defined clouds corresponding to the trajectories in the abovementioned clusters (Fig. 9a). The fast Antarctic transport cluster 2 shows much lower ice formation efficiency. While the Australian cluster 1 shows highest ice formation efficiency, except at the temperature range −25 to −20 ^∘C. There, the marine clusters 3 and 4 exhibit larger ice formation efficiencies, albeit within statistical uncertainty.

The MACC model was used to obtain simulated aerosol conditions for each classified well-defined cloud (see Sect. 2.3). The aerosol effect on ice formation efficiency in the observed clouds is investigated by categorizing these clouds based on the simulated aerosol mixing ratios. As discussed in Sect. 3.2, the dominance of sea salt in the simulated aerosol mixing ratios above Lauder might obstruct the assessment of aerosol impact on the formation of mixed-phase clouds. Therefore, two different quantities were defined based on the MACC aerosol mixing ratios to use as a proxy for aerosol abundance. The first one is the MACC total aerosol mixing ratio excluding the sea salt mixing ratio (see the time series in Fig. 5). The second one is the former divided by the sea salt mixing ratio, i.e., the MACC aerosol mixing ratio excluding sea salt relative to sea salt. Varying arbitrary thresholds on these two quantities were then used to categorize the abovementioned clouds. Figure 10 shows the fraction of ice-containing clouds against cloud-top temperature for all well-defined clouds and for well-defined clouds corresponding to increasing MACC total aerosol mixing ratios excluding sea salt as well as increasing MACC aerosol mixing ratios excluding sea salt relative to sea salt. Enhanced ice formation efficiencies were found for clouds experiencing a combination of high aerosol mixing ratios (excluding sea salt) and a low sea salt fraction. This is an indication for higher ice formation efficiency for clouds which are more strongly influenced by continental aerosol.

Figure 10 Same as Fig. 8 but for all well-defined clouds (black, same as in Fig. 8) and for well-defined clouds corresponding to a combination of criteria, namely, MACC total aerosol mixing ratio excluding sea salt relative to sea salt > 0.3 and MACC total aerosol mixing ratio excluding sea salt > 0.3 µg kg⁻¹ (orange), as well as > 0.5 and > 0.6 µg kg⁻¹, respectively (violet).

4 Summary, conclusion, and outlook

This study investigated the impact of variations in the atmospheric aerosol load on the frequency of heterogeneous ice formation over Lauder, New Zealand / Aotearoa. As prerequisites, long-term datasets of polarization lidar observations, atmospheric temperature and humidity profiles from both radiosondes and model reanalyses, and HYSPLIT backward trajectories, as well as of CAMS-MACC aerosol reanalysis data, were acquired and processed. An 11-year polarization lidar dataset from Lauder was analyzed with the Seifert et al. (2010, 2015) methods to assess the ice formation efficiency in clouds with top temperatures in the heterogeneous freezing range from −40 to 0 ^∘C. Thereby, the cloud dataset considered only clouds with a lidar-detectable cloud top height, in order to prevent positive temperature biases by deep, attenuating clouds with underestimated cloud tops.

Based on a long-term radiosonde dataset from Invercargill / Waihōpai, New Zealand / Aotearoa, the general conditions for the formation of mixed-phase clouds in the region of southern New Zealand / Aotearoa were characterized. It was found that the −35 ^∘C isotherm lies between 7 km (austral winter) and 8.4 km (austral summer). In the heterogeneous freezing temperature range, dew point spreads of 0–2 K occurred with a frequency of 14 % (austral summer) and 17 % (austral winter).

Next, the general aerosol conditions over Lauder were assessed based on a statistical analysis of HYSPLIT back trajectories with the TRACE tool, as well as based on a long-term statistics of CAMS-MACC aerosol reanalyses data. Generally, it was found that sea salt and marine air masses are the dominating aerosol species. Nevertheless, this constant marine background was superimposed by aerosol of considerable variability regarding the air mass source regions and the aerosol load. Both clean, pristine air masses from Antarctica and the Southern Ocean and aerosol-laden air masses from continental Australia occur frequently over New Zealand / Aotearoa. The number of identified cases per air mass was sufficient to split the lidar-based cloud dataset into different subpopulations in order to assess the impact of air mass origin on the frequency of heterogeneous ice formation.

Overall, ice formation efficiency in clouds above Lauder was found to be lower than at Leipzig in the polluted midlatitudes of the Northern Hemisphere but higher than in very clean and pristine environments, such as Punta Arenas in southern Chile, 8000 km downwind of any neighboring landmass. This can be expected as New Zealand / Aotearoa is influenced both by Australian and marine aerosol. Further, the incorporation of the trajectory tools and aerosol model simulations both point clearly to higher ice formation efficiency for clouds which are more strongly influenced by continental aerosol and to lower ice formation efficiency for clouds which are more influenced by Antarctic/marine aerosol and air masses.

In conclusion, it can be stated, based on the outcome of the study, that the atmospheric conditions over the South Island / Te Waipounamu of New Zealand support contrast studies for evaluation of the relationship between aerosol load and mixed-phase cloud formation. We argue that this region is one of the more accessible places in the Southern Hemisphere midlatitudes, where local aerosol sources are scarce while the general circulation permits alternating periods of pristine and aerosol-laden air masses, respectively.

Our study made use of a unique 11-year polarization lidar dataset from Lauder. Other active remote sensing instrumentation was not available for such a long time period over this region. There is, hence, a high potential for future long-term experiments that incorporate additional observational capabilities, specifically of Raman and Doppler lidar as well as cloud Doppler radar. Radenz et al. (2021a) have extensively demonstrated the capabilities of these additional instruments in characterizing (i) the ice-nucleating capabilities of the aerosol mix at cloud level, (ii) the impact of atmospheric gravity waves on the sustained formation of supercooled liquid water, and (iii) the determination of aerosol impacts on ice water mass production as estimated from cloud radar observations. Moreover, cloud radar would allow us to observe the cloud layers in their full extent, which would minimize the number of clouds with undefined cloud top heights/temperatures.

Future studies might also make further use of other atmospheric measurements from the Lauder observatory. Profiles of trace gas concentration could be used as indicators for Southern Ocean and marine air masses. From February 2023, the aerosol lidar at Lauder has been automated for 24 h operation, measuring in alternating 10 min intervals. This will greatly increase the data available for possible extension or refinement of this study.