GEOS-GOCART simulation of IMO impact on aerosol fields
All simulation experiments were run with the Goddard Chemistry Aerosol Radiation and Transport (GOCART) aerosol module39,40 in NASA Goddard Earth Observing System (GEOS) Earth System Model (ESM). The GEOS model has a one-moment cloud microphysics module and a rapid radiation transfer model for general circulation models (RRTMG). Sea surface temperature (SST) for the atmospheric dynamic circulation is provided by the GEOS Atmospheric Data Assimilation System (ADAS) that incorporates satellite and in situ SST observations. The model is run in the replay mode using meteorological fields from the Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) reanalysis41. “Replay” mode sets the model dynamic state (winds, pressure, and temperature) every 6 h to the balanced states provided by the Meteorological Reanalysis Field of the Modern Research and Applications Reanalysis Version 2 (MERRA-2). We run GEOS at a global horizontal resolution of approximately 50 km on a cubic sphere grid and 72 vertical layers from the surface to 0.01 hPa. The time step for dynamic calculation is 450 s. The temporal resolution of the radiation is 1 h. All experiments run from 201910 to 202012, with the first three months as the spin up period. We use monthly results in our estimation of forcing.
We have two set of experiments: business as usual (BAU) and Covid impact (Covid) emissions. The BAU used anthropogenic emissions of aerosols and precursor gases from the Community Emission Data System (CEDS)42 but repeat the 2019 emissions for 2020. The dataset includes nine emission sectors (energy, industry, road transportation, residential, waste, agriculture, solvent, shipping, and air traffic). Biomass burning emissions were taken from the GSFC-developed Quick Fire Emission Dataset (QFED)43. Volcanic emissions come from the dataset that is based on the satellite volcanic SO2 observations from the OMI instrument on board the Aura satellite. Biogenic emissions were calculated with the Model of Emissions of Gases and Aerosols from Nature (MEGAN) that is embedded in GEOS model. The wind-driven emissions, such as dust and sea salt, were calculated on-line. Time varying greenhouse gases, such as CO2, CH4, N2O and ozone-depleting substances, were provided by CMIP5 project.
The second set of experiments (Covid) adjusted BAU anthropogenic emission to reflect the impact of Covid. Using mobility data from Google and Apple44, daily scale factors in 2020 were derived on sector bases for ten species including important aerosols and their precursors. Because of the rapidly changing emissions due to various timing and strength of lockdown measures, daily scale factors were provided not only to scale down emission amounts but also to move emissions from monthly to daily. The Covid adjusted daily anthropogenic emissions were generated by applying these scale factors to CEDS 2019 monthly emission.
A summary of S emissions under different scenarios is provided in Table S1. For each set of experiments, there are three scenarios: full ship emissions, reduced ship emissions using the IMO2020 standards, and no ship emission of S. Other emissions are kept the same. We take the difference in aerosol loading between with full ship emissions of SO2 and with reduction due to IMO2020 as the impact of IMO2020. The difference in aerosol loading is translated into Nd changes with method in the following subsection.
Deep learning-based Nd
The operational version of NASA’s Global Earth Observing System (GEOS) runs single moment cloud and aerosol microphysical schemes. They do not predict cloud condensation nuclei (CCN) and cloud droplet number concentrations (Nd). We estimate Nd using a diagnostic deep learning-based approach, involving the usage of two neural network (NN) parameterizations. The first NN (termed MAMnet) is an emulator for the Modal Aerosol Module, which takes bulk aerosol mass for 5 externally-mixed species (sulfates, sea salt, dust, black carbon, and organics) and the atmospheric state (temperature, pressure) as input, and predicts the number concentration and composition for 7 internally-mixed lognormal modes (accumulation, aitken, coarse/fine dust, coarse/fine sea salt, primary carbon matter). MAMnet was trained on 5 years of data from a GEOS simulation implementing the MAM7 aerosol module, and validated against ground observations. The CCN concentration at a given supersatureation can be readily estimated using the 7-modal size distribution and composition45. A second NN model (Wnet) is used to estimate Nd. Estimation of Nd requires the characteristic vertical wind velocity (W) at the scale of individual parcels, typically proportional to its subgrid-scale standard deviation (σW). Wnet takes the atmospheric state, as well as coarse metrics of turbulence (Richardson number and total scalar diffusivity) as inputs and predicts σW for each grid cell. Wnet was trained on 2 years of a global, non-hydrostatic, storm-resolving simulation of the GEOS model and physically constrained by ground-based observations of σW from around the world using a novel generative approach46. The aerosol size distribution predicted from MAMnet as well as σW are used to predict Nd using the Abdul-Razzak and Ghan scheme47. Our NN-based method emphasizes observations and conservation of known physics during the development of the NNs and ensures a robust prediction of CCN and Nd.
Calculating aerosol indirect forcing
We use the same methodology reported in Yuan et al.21. We consider the Twomey effect and effects of cloud fraction and LWP adjustments. Without considering aerosol effects on cloud fractions, cloud albedo sensitivity to aerosols can be taken as the sum of the Twomey effect and aerosol induced LWP adjustments:
$$S=\fracdA_cdN_d=\fracA_c(1-A_c)3N_d\times \left(1+\frac52\fracdlnLWPdlnN_d\right)$$
(1)
where S is the susceptibility of cloud albedo (Ac) to droplet number concentration \(N_d\)16.
We then have
$$\varDelta SW_TOA=-SW_downwelling\times Cf\times {{\rmS}}\times \varDelta N_d$$
(2)
Aerosol indirect forcing from the Twomey effect and LWP adjustment is therefore:
$$\varDelta SW_TOA= -SW_downwelling\times Cf\times A_c\times (1-A_c) \\ \times \left(\frac13+\frac56\fracd\mathrmlnLWPdlnN_d\right)\times \varDelta \mathrmlnN_d$$
(3)
To consider the effect of Cf adjustment due to aerosols, we consider the sensitivity of scene albedo (A) to Nd. \(A=A_acCf_total+A_s(1-Cf_total)\). We have:
$$S^* = \fracdAdN_d=\fracd(A_acCf_total+A_s(1-Cf_total))dN_d \, \approx \, Cf\times S \\ +\left(1-Cf_high\right)\times \fracdCfdN_d\times (A_ac-A_s)$$
(4)
where A is the scene albedo, i.e., including both cloudy, Aac, and clear, As, parts; \(Cf_total\) and \(Cf\) are all cloud and low cloud fraction obtained from the MYD08_M3 data; As is the surface albedo, derived from the CERES EBAF-TOA data48; \(1-Cf_high\) is used to take into account of effect of overlap on Cf adjustment. We assume a maximum overlap between high and low clouds. We assume minimum aerosol effects on high clouds. The estimation is done at monthly time scales.
CF and LWP adjustments, \(\fracdCfdN_d\) and \(\frac{d{{{\mathrmln}}}LWP}dlnN_d\), are derived from our previous work based on large number of ship-track samples21. The assumption is that clouds with similar properties respond similarly to addition of aerosols and ship-tracks detected under diverse background cloud conditions can be used to effectively derive these adjustments. Our results are based on the responses from observed ship-track sampled under diverse set of environmental conditions, which allows us to derive robust cloud adjustments based on numerous ship-track samples21. There are a few assumptions and approximations as noted in our previous study21 and we reiterate them here. We used SWdownwelling at the surface from CERES instead of SWdownwelling at the cloud top, which underestimates the total forcing since SWdownwelling at the cloud top is larger. The LWP and Cf adjustments can be sensitive to more variables that those considered here. We assume the derived \(\fracdCfdN_d\) and \(\frac{d{{{{{\mathrmln}}}}}LWP}dlnN_d\) and their dependence on background variables apply to regions that have less ship-track samples. Also, potential semi-direct effects due to absorbing aerosols from ship-emissions are not explicitly addressed in this study since we do not have enough observations.
The cloud adjustments used here can depend on the background cloud Nd, SST, EIS, and background Nd and thus they have spatiotemporal variations due to background changes. The dependence of cloud adjustments can also be parameterized with one or more variables. The Nd-only functional form provides a lower bound on the forcing calculation21 and we explore different 2-variable combinations to provide a range of estimates. We combine the cloud adjustments with the simulated \(\varDelta N_d\) due to IMO2020 and observations of clouds and other parameters in 2020 to calculate its forcing. We report the mean forcing from all five functional forms of cloud adjustments as well as the standard deviation. We also report the warming effect expected from both the upper and lower bounds in Fig. 3. Due to the Nd-dependent nature of cloud adjustments, the same \(\varDelta N_d\) can result in different magnitude of radiative forcing21,49.
There is systematic difference between GEOS-modeled and MODIS observed climatology of Nd. At each grid point, we calculate the ratio between modeled and observed Nd based on monthly data and scale the modeled \(\varDelta N_d\) with the ratio before using above equations to calculate the forcing. The global mean values change by 10% between scaled and non-scaled \(\varDelta N_d\), all coming from the CF adjustment since the LWP adjustment and the Twomey effect depend on \(\varDelta N_d/N_d\) that does not change with scaling. Regionally, the difference can be as large as 30%, e.g., in the Southeast Atlantic.
The CERES EBAF-TOA data48 provides monthly and climatological averages of observed top-of-atmosphere and computed cloud radiative effect and absorbed solar radiation. The top-of-atmosphere net fluxes provides constraints to the ocean heat storage. It is used here to calculate the interhemispheric contrast in absorbed solar radiation and energy balance. We note that although the interhemispheric contrast is a residue of two large numbers, e.g., the amount of mean absorbed solar radiation in both hemispheres, the observed variation of the contrast is always small. Therefore, even though we cannot directly attribute the variations in the interhemispheric contrast to IMO 2020, it is reasonable to discuss their temporal evolutions and compare the IMO 2020 impact with the observed changes.
Transient warming of IMO2020
We consider the simple one-layer energy balance model27:
$$C\times \fracdTdt=F-\lambda \times T$$
(5)
where C is the heat capacity of the well-mixed ocean layer, T is the temperature anomaly from the equilibrium, t is time, F is the forcing, and \(\lambda\) is the climate feedback parameter. For an abrupt forcing, the solution is:
$$T=(F/\lambda )\times (1-e^-\lambda t/C)$$
(6)
Using C = 8.2 W yr/m2/K and \(\lambda\) = 1.2 Wm−2K−1, for a forcing of F = 0.2Wm−2, we get the temperature change at the new equilibrium is 0.17 K with a time scale of C/ \(\lambda\) = 7 years. The warming rate is F/C = 0.2/8.2 K/yr = 0.024 K/yr or 0.24 K/decade. \(\lambda\) has uncertainty associated with it and its 1-\(\sigma\) is 0.25 Wm−2K−1. With this, we can estimate equilibrium \(T\) to be between 0.14 and 0.21 K. Equation 6 is used to calculate the expected warming trajectory in the 2020 s when combined with a simple long-term upward trend in Fig. 3. The observed global mean temperature is from the National Aeronautics and Space Administration (NASA) Goddard Institute of Space Studies.
Contributions of background Nd, CF, and \(\varDelta N_d\)
We calculate the annual mean of incoming solar radiation for each oceanic grid in the North Atlantic and use this map of seasonally invariant incoming solar radiance to calculate IMO2020 forcing (see section c of Methods). The seasonal cycle of forcing using the seasonally invariant solar radiation is shown in Figure S1, which differs substantially from Fig. 3a, highlighting the impact of seasonal cycle in solar radiation. The peak season for the forcing is now wintertime instead of summertime. This serves as our baseline to test sensitivity of forcing to different variables.
The sensitivity of the forcing to each factor is assessed through the following procedure. We first calculate the seasonal variations of the forcing using observations that contain its seasonal variations. We then calculate a map of annual mean for each variable and use it to calculate the forcing, effectively removing its impact on the seasonal changes. The relative difference between these two calculations can be taken as a measure of how much each variable contributes to the seasonal changes.
