The impacts of shipping emissions on lightning: roles of aerosol-radiation-interactions and aerosol-cloud-interactions

We apply global observations of lightning strokes from the ground-based World Wide Lightning Location Network (WWLLN) (Virts et al 2013). WWLLN locates lightning strokes using radio receivers distributed across the world. The observation network was initiated in 2003 with 11 sensors (Lay 2004), with the number of sensors increasing to more than 70 in 2013, providing the world's longest record for global lightning observations. These sensors detect impulsive signals from lightning discharges that have steadily propagated in the Earth–ionosphere waveguide within a time window of one microsecond, and each sensor can cover lightning strokes within a distance of about 3000 km (Dowden et al 2002). WWLLN detects very low-frequency ratio (VLF) in lightning and preferentially detects cloud-to-ground (CG) strokes due to its higher VLF value (Rodger et al 2009, Virts et al 2013). Comparison of WWLLN observations with satellite-based observation systems (e.g. the lightning imaging sensor and the optical transient detector satellite instruments, LIS/OTD) and local lightning detection network at overlapping observational years suggests that these observational networks are generally consistent in the magnitude and spatiotemporal distribution of lightning frequencies (Lay 2004, Jacobson et al 2006, Bürgesser et al 2017).

We follow Thornton et al (2017) to focus on lightning activities over the northeastern Indian Ocean, which covers one of the most heavily trafficked shipping lanes in the world. Figure 1 shows the shipping emissions of fine particulate matter (PM_2.5) from the Emissions Database for Global Atmospheric Research (EDGAR) and the observed annual mean lightning frequency over this region averaged over 2005–2021. The WWLLN observation clearly shows the hotline of lightning frequency along the shipping lane region centered at around 6° N (5–7° N, 82–92° E), compared to an adjacent reference region (7–9° N, 82–92° E) covered by comparable amounts of sensors (figure 1(b)). Figure 1(c) further shows that the lightning frequency over shipping lanes is about 20%–100% higher compared with that in the reference region during 2005–2021. Such difference is especially notable in years 2013, 2015, and 2018. Since the two oceanic regions are geographically close to each other, we propose that the large difference in lightning frequency between the two regions is more likely due to shipping-induced aerosols which modulate the strength of convection, rather than large-scale meteorological conditions. We also see an increasing trend of lightning activities in the shipping lane region that might be driven by increasing ship emissions, as also proposed by Thornton et al (2017), but this may partly reflect the trend in the detection capability with the increasing number of sensors.

Zoom In Zoom Out Reset image size

We use the Weather Research and Forecasting model coupled with Chemistry (WRF‐Chem) version 3.6.1 to interpret lightning activities and the impact of shipping emissions. WRF-Chem is an online-coupled meteorology-chemistry model (Grell et al 2005, Fast et al 2006, Powers et al 2017), in which aerosols released and/or chemically generated from shipping emissions can interact with radiation as well as the cloud thermodynamics and microphysics. We design three nested modeling domains centered on the shipping lanes in the northeastern Indian Ocean as shown in figure 1(a). The horizontal resolutions and grid numbers for the three domains are 36 km with 112 × 111 grids, 12 km with 277 × 118 grids, and 4 km with 394 × 172 grids, from the outermost and innermost domains, respectively. We set 50 vertical layers from surface to top pressure at 50 hPa for each domain. Meteorological initial and boundary conditions are provided by the National Centers for Environmental Prediction FNL (Final) Operational Global Analysis data at 1° × 1° resolution. We acknowledge that using the FNL data as initial meteorological fields has inevitably introduced some impacts of shipping emissions on initial background meteorological conditions by assimilating in-situ and remote sensing observations. However, as our study focuses on the short-term responses of convections and lightning activities to shipping emissions, we suggest that the initial meteorological conditions should not significantly affect such response if the model spin-up time is sufficiently long compared to the lifetime of a typical concentration system. We also do not apply any nudging techniques in the simulation to avoid the influence of data assimilation on the results during the simulation period. The initial and lateral boundary chemical conditions are archived from the global simulation using the model for ozone and related chemical tracers model.

The physical and chemical schemes used in our WRF-Chem simulation are listed in table S1. Specifically, we use the 8-bin version of the model for simulating aerosol interactions and chemistry (MOSAIC) for aerosols including sea salts and dust (Zaveri et al 2008), the Lin scheme for microphysics (Lin et al 1983). The MOSAIC module assumes aerosols to be internally mixed so that all particles within a bin have the same chemical composition, while particles in different bins are externally mixed (Zaveri et al 2008). We use the Grell 3D ensemble cumulus scheme (Grell and Devenyi 2002) for the outermost two domains with resolution coarser than 10 km where a convective parameterization is necessary, and do not use the cumulus scheme for the innermost domain since the convective processes can be explicitly solved at a resolution of 4 km (Kain et al 2008). The parameterization of lightning activity is closely related to the convective process that will be described in section 2.3.

Anthropogenic emissions are from the EDGAR (https://edgar.jrc.ec.europa.eu/) inventory. EDGAR provides monthly mean shipping emissions from 1970–2018 at 0.1° × 0.1° spatial resolution for SO₂, NO_x , CO, total non-methane volatile organic compounds (NMVOC_s), NH₃, as well as particulate matter compounds, including PM_2.5, PM₁₀, black carbon (BC), and organic carbon (OC) (Crippa et al 2016). The shipping emissions in EDGAR are derived based on real-time shipping information from the automatic identification system (AIS) and emission factors (Lack et al 2009). The emissions inventories are then a sum over emission factors scaled by the fuel usage statistics derived from the shipping activity databases. Nevertheless, uncertainties in the shipping information (tracks, fuel type, vessel type, etc) recorded in the AIS database and emission factors may contribute to the uncertainties in the EDGAR shipping emission inventory (Russo et al 2018, Mao et al 2020). The EDGAR inventory also does not include the daily variation of shipping emissions. This limits the model's ability to capture the daily variations in shipping's impact on lightning activity. The model does not implement an online plume model for the shipping emissions (Mao et al 2021).

Figure 1(c) shows that 2013, 2015, and 2018 are the years with the highest lightning frequency over the shipping lane of our interest. Examination of its monthly variation suggests that the lightning activities in this region peak in the boreal winter due to the southward mitigation of the intertropical convergence zone (figure S1). The shipping emissions in January 2013, 2015, and 2018 are 92.9%, 99.9%, and 107.3% compared to the 2010–2018 mean, suggesting that the high lightning frequency in these years is not solely due to differences in shipping emissions. We thus choose three periods with relatively high lightning flash counts scattered in January 2013 (4–6 January), 2015 (25–29 January), and 2018 (28–31 January) as the cases for WRF-Chem simulation (figure 1(d)). We follow previous studies (Federico et al 2014, Gharaylou et al 2020, Guo et al 2022) to set the first 12 h in each simulation for spin-up to allow the meteorological conditions to be adjusted to different scenarios.

The lightning parameterization used in WRF-Chem is based on cloud-top height ($z$) (Price and Rind 1992, hereafter PR92), which has been widely implemented in atmospheric chemistry models (Wang et al 2011, Murray et al 2012). Based on the framework of Vonnegut (1963) and Williams (1985), the lightning frequency in PR92 parameterization is established to be proportional to the fifth-power of convective cloud-top height. Michalon et al (1999) further modified the equation from PR92 for marine flash (${f_m}$):

$$\begin{equation}{f_c} = 3.44 \times {10^{ - 5}}{z^{4.9}}\end{equation} \tag{ 1 }$$

$$\begin{equation}{f_m} = 6.57 \times {10^{ - 6}}{z^{4.9}}.\end{equation} \tag{ 2 }$$

The WRF-Chem model determines the cloud-top height $z$ based on the level of neutral buoyancy, which is calculated based on the Grell–Devenyi convective parameterization. In the innermost domain where the cumulus scheme is turned off, we determine the top level with radar reflectivity larger than 20 dBZ as $z$ instead. A fixed but region-varied intra-cloud (IC):CG ratio is applied to separate IC and CG flashes in the WRF-Chem model according to Boccippio et al (2001) to match lightning measurements.

The model uses equations (1) and (2) to calculate the flash rate for continental and marine clouds in every timestep ($\Delta t$), and returns the accumulated lightning counts ($A$) in each model grid as

$$\begin{equation}{A_{t + 1}} = {A_t} + f\Delta t.\end{equation} \tag{ 3 }$$

However, since the model calculates total lightning counts instead of lightning density in each model grid, the simulated frequency derived from lightning counts is grid-scale dependent, i.e. the lightning frequency is bound to the resolution (area) of the model grid. We follow Wong et al (2013) to apply a factor ($c$) introduced by Price and Rind (1994) to the simulated flash rate at each model grid. The factor is a function of model grid size, given as:

$$\begin{equation}c = 0.97241{e^{0.048203R}}\end{equation} \tag{ 4 }$$

where $R$ is the grid area in squared degrees, and $c$ represents the resulting factor scale to the flash rate at the given grid resolution. We note that, however, model at grid resolution within 4–6 km is less influenced by this grid-scale dependence since the initial cloud-top height parameterization is based on radar data at a resolution of 5 km (Williams 1985).

Our Base simulation applies the above settings that include shipping emissions and fully considers both ARIs and ACIs. We also conduct seven sensitivity simulations to assess the influence of shipping emissions on lightning frequency separately through ARIs and ACIs, by excluding shipping emissions in the model simulation and inactivating ARIs and/or ACIs (table S2).

Simulations A (Base simulation) and B, with both ARIs and ACIs available but shipping emissions turned on and off, respectively, quantify the integrated influence of shipping emissions on lightning frequency. ACIs are simulated by linking prognostic aerosols to cloud condensation nuclei (CCN). We use a prescribed cloud droplet number concentration of 10 cm⁻³ to exclude ACIs in sensitivity simulations C and D, following previous studies (Zhao et al 2017, Zhang et al 2018, Zhou et al 2019). This intends to assess the contribution of ACIs to lightning frequency from shipping emissions. ARIs are simulated by linking aerosol optical properties (such as aerosol optical depth (AOD), single scattering albedo, and asymmetry factor) to the radiation scheme (Chapman et al 2009). We conduct simulations with ARIs turned off in the presence (simulation E) and absence (simulation F) of shipping emissions to investigate the contribution of ARIs to lightning frequency from shipping emissions. Absorbing aerosols (such as BC) also influence cloud distribution by modifying atmospheric temperature profiles, known as the aerosol semi-direct effect. In this study, ARI's impact quantified as the difference between simulations E and F also included the aerosol semi-direct effect (Zhang et al 2018). We additionally conduct simulations G1/G2 by inactivating ARIs and ACIs in the presence/absence of shipping emissions, in order to verify whether shipping emissions will influence lightning frequency if there is no aerosol-climate feedback. All these simulations are conducted for the three selected cases in 2013, 2015, and 2018.

Figures 2 and S2–S4 evaluate the simulated AOD, daily precipitation, and lightning flash counts in the shipping lanes region (5–7° N, 82–92° E) in the WRF-Chem model for the three cases. The datasets for evaluation are the second version of Modern-Era Retrospective Analysis for Research and Applications for AOD, the Tropical Rainfall Measurement Mission (TRMM) satellite for precipitation, and WWLLN for lightning flash rate, respectively.

Zoom In Zoom Out Reset image size

Figures S2–S4 show shat our WRF-Chem simulation captures the overall spatial distributions of AOD and precipitation over the region of interest for all three cases. with the spatial correlation coefficients of 0.64 (N= 12 210, p< 0.01) and 0.44 (N= 32 292, p< 0.01) in the 2015 case as an example for AOD and precipitation, respectively. It reproduces the high AOD values over the Indian Peninsula and Indochina Peninsula, but overestimates that over the Bay of Bengal (figure S2). The model also reproduces the high precipitation with over 100 mm per day over Indonesia, consistent with the TRMM observations, but overestimates precipitation over the oceanic region at around 6° N, 90° E.

Capturing the magnitude and temporal evolution of lighting flash counts in meteorology-chemistry coupled models is much more challenging than that of AOD and precipitation, as it involves the combined uncertainty in lightning parameterization and in simulating convection, aerosol properties and its interaction with radiation and clouds. Our WRF-Chem simulation largely reproduces the temporal variation of the WWLLN observed lightning counts for all three cases (figures 2(a)–(c)), indicating the simulation successfully reproduces the temporal evolution of convections and lightning activities. However, it overestimates lightning count by an order of one or two magnitudes. Similarly, Wong et al (2013) also found that the WRF simulation biased high in lightning flash rates by an order of magnitude compared to ground-based detection networks in the United States. This discrepancy reflects both measurement bias and model deficiency. WWLLN sensors preferably detect CG flash (section 2.1), but have limited ability to detect IC flash. Even though the current WWLLN network has been shown to capture more lightning counts than the satellite-based LIS/OTD dataset in the equatorial Indian Ocean (Virts et al 2013), it may still underestimate the true frequency of lightning counts due to an insufficient number of sensors and its inability to count IC flashes. From the model perspective, the PR92 scheme solely relies on the cloud top height as an indicator of convention strength, however, stronger convections are not always indicative of more frequent lighting activities (Xu et al 2017). The high bias in the simulated lightning counts may also reflect the model's overestimation of precipitation over the shipping lane region (figures S2–S4). More efforts are required to improve the accuracy of the simulation of lightning in weather or climate models.

We find that lightning counts averaged over the shipping lanes (82–92° E, 5–7° N) in the Base simulation (simulation A) are 2.7%, 5.5%, and 4.4% larger than those in the simulation without shipping emissions (simulation B), in the simulation periods of 2013, 2015, and 2018, respectively. Such difference is particularly large when convection is active and lightning frequency is high, reaching 22.3% on 28 January 2015 and 10.5% on 29 January 2018 (table 1). We also note that the differences in lightning frequency introduced by shipping emissions show temporal variabilities, with larger lightning frequency in some periods when shipping emissions are excluded (figures 2(d)–(f)).

The shipping-induced enhancement in the lightning counts can be entirely attributed to the interactions between shipping-induced aerosols and the atmosphere through ARIs and ACIs. This is supported by our simulations of G1 and G2 (table S2), which show that the simulated lightning counts remain exactly the same when shipping emissions are included but with both ARIs and ACIs invalidated in the simulation. We quantify that shipping activities enhance the AODs by 1.3%–7.2% over the shipping lane region in the three cases by direct emissions of particles (PM_2.5, PM₁₀, BC, OC) and secondary formation of aerosols from the emission of precursors (SO₂, NO_x ). The shipping-induced enhancement of aerosols, even at a relatively small magnitude, has shown moderate influence on the mean lightning frequency in the three cases.

We next examine the involved physical processes of how these shipping-induced aerosols alter the lightning frequency. We analyze the regional mean difference in simulated cloud droplet number, cloud water content, radar reflectivity in the altitude range of 0–20 km, and cloud-top height attributable to shipping emissions during the period with the most active convection in each case. We find that shipping-induced aerosols increase the convection strength and the lightning frequency by 6.2%–22.3% for all three cases, when ACIs and ARIs are both considered (table 1). Mechanically, the shipping-induced aerosols serve as additional CCN and increase the cloud droplet number by 3.2%–8.9%, leading to more but smaller cloud droplets according to the Twomey effect (Twomey 1977). With the reduction of the effective radius of cloud droplets, the collision-coalescence processes are hindered and the precipitation of cloud water is delayed. Our simulations suggest both an invigoration of updrafts (0.3%–4.8% increase in cloud top heights) and cloud water content (1.2%–6.5%) (table 1, figure 4). That is to say, the enhancement in lightning, although parameterized here using only cloud top height, may be both dynamic and microphysical in nature. While liquid and ice water content are altered by the Twomey effect described above, the invigoration of updrafts may result from some combination of (1) increased latent heat release while freezing (Rosenfeld et al 2008) and (2) increases in environmental humidity (Abbot and Cronin 2021). The latter is indicated by small increases in environmental humidity due to shipping-induced aerosols in all three cases (figure S5).

The increase in regional mean cloud-top height of 0.3%–4.8% in response to the consistent increase in cloud water and radar reflectivity is responsible for the enhancement of lightning frequency, however, it is not proportional to the increase in lightning flash count. This may reflect the nonlinear relationship between cloud-top height and flash count in the PR92 parameterization. Our analyses thus support that shipping emission aerosols foster the convection through microphysical processes and lift the convective cloud-top height, ultimately contributing to the enhancement of lightning counts.

The microphysical processes and their contribution to lightning activity also involve strong temporal evolution, as illustrated in the 2015 case (figure 3). The most significant convection events occurred at around 0 UTC and 22 UTC on 27 January (figure 3(e)). We find that shipping emissions contribute to substantial increases in the mixing ratio of the cloud droplet number, cloud water content, and radar reflectivity in these two periods, effectively enhancing the lightning frequency. In comparison, there is an apparent decrease in cloud water content and radar reflectivity before the activation of convection. This is because when additional aerosols are activated as CCN, the effective radius of cloud droplets becomes smaller, which suppresses the formation of precipitation and convection at its early stage, according to the Twomey effect. Cases in 2013 and 2018 also show similar patterns (figures S6–7). This also explains the temporal variability of shipping-induced lightning flash counts as shown in figure 3.

Zoom In Zoom Out Reset image size

We use the methods described in section 2.4 to separate the influence of ACIs and ARIs on lightning frequency, and results are summarized in table 1. We find that when ACIs are inactivated in the model, shipping emissions lead to changes in lightning flash counts of +4.5%, +0.25%, −6.6% for the 2013, 2015, and 2018 cases, respectively. These responses are significantly lower than the consistent increase of 6.2%–22.3% when both ACIs and ARIs are considered. Such differences are also reflected in the simulated cloud droplet number, cloud water content, radar reflectivity, and cloud-top height. Since the cloud droplet has been set to 10 cm⁻³ in the model to inactivate ACIs, the simulated CCN is fixed to a relatively low concentration, limiting the role of the shipping-induced aerosols to serve as CCN and to develop the convective cloud via the microphysical processes as described in section 3.2. This is supported by a much smaller shipping-induced increase in the radar reflectivity of only about 0.1% when the ACIs are inactivated, compared to the 0.2%–2.1% increase when both ACIs and ARIs are considered in the model, and the opposite contribution to radar reflectivity during the active convection period when excluding ACIs in all cases as shown in figure 4. In response, the cloud top height decreases by −2.3% and −1.2% in 2015 and 2018 cases with inactivated ACIs. The above results support the importance of ACIs in invigorating convections and enhancing lightning frequency.

Zoom In Zoom Out Reset image size

In comparison, ARIs' impacts on lightning frequency are inconsistent in the three cases. In the 2018 case, the removal of ARIs increases the shipping-induced lightning enhancement from 10.5% to 15.8% (table 1), indicating that ARIs suppress convection and lightning activity, which is consistent with Zhang et al (2020). In contrast, ARIs enhance convection and lightning activity in 2013 and 2015 cases, and the contributions are comparable to ACIs. We also notice that inactivating either ACIs or ARIs both show much fewer changes in lightning counts (+0.25% and −2.6%) compared to the 22.3% increase when both ACIs and ARIs are active in the 2015 case, suggesting that there can be strong interactions between ACIs and ARIs (Lohmann and Feichter 2001). Figure 4 further shows inconsistent or uncertain responses of cloud properties and lighting activities to shipping emissions when ARIs are excluded from the model.

We find the inconsistent responses of lighting activities to ARIs are associated with different aerosol levels in three selected cases. Aerosol loading over the shipping lanes was relatively low in 2018, as evidenced by a 50% lower regional mean AOD in 2018 than in 2013 and 2015 cases. With higher aerosol loadings and with ARIs activated in the simulation, shipping-induced aerosols heat the air temperature at around 3–5 km and cool the air temperature at around 1 km above the surface in the 2013 and 2015 cases (figure S8). This may reflect the radiative effect of shipping-induced sulfate and BC aerosols. Sulfate aerosols are accumulated below 1 km above the surface (figure S9). These aerosols could scatter the shortwave radiation to the atmosphere ahead (3–5 km above the surface), where BC aerosols there further absorb the radiation and heat the atmosphere. Moreover, ingestion of aerosols in warm clouds (which can reach around 5 km according to figure 3(c)) may increase the albedo and the cloud reflectance (Twomey 1977), resulting in the heated air temperature at around 3–5 km. The anomaly heating then leads to a positive anomaly in the vertical velocity at 3–5 km (red dashed line in figures S8(a), (b), and hence invigorates convection aloft and enhances lightning flash counts.

We hypothesize that the low aerosol loadings in the 2018 case may limit the scattering effect and the resulting heating of the lower atmosphere as in the 2013 and 2015 cases, and hinder the development of convection as supported by the much smaller change of vertical velocity in Figure S8c compared to figures S8(a) and (b). Furthermore, we find that the average cloud droplet number over the shipping lanes in 2018 is much lower than that in 2013 and 2015 (4.7 × 10⁵ kg⁻¹ versus 1.1–1.4 × 10⁶ kg⁻¹). It is possible that the heating and evaporation of cloud droplets may significantly reduce the cloud droplet numbers, impeding the ACIs-involved cloud microphysics and the development of convection (Liu et al 2022). The evaporation itself can also suppress the convection because of reduced latent heat. The ARIs impact also includes the influence on cloud status and distributions through altering the atmospheric temperature profiles, but which is difficult to independently quantify. These results point to a complicated effect of ARIs on convection that might be tied to the number, composition, and vertical structure of aerosols, as well as the interaction with ACIs.