Substantial Cloud Brightening From Shipping in Subtropical Low Clouds

The influence of aerosol particles on cloud reflectivity remains one of the largest sources of uncertainty in our understanding of anthropogenic climate change. Commercial shipping constitutes a large and concentrated aerosol perturbation in a meteorological regime where clouds have a disproportionally large effect on climate. Yet, to date, studies have been unable to detect climatologically relevant cloud radiative effects from shipping, despite models indicating that the cloud response should produce a sizable negative radiative forcing (perturbation to Earth's energy balance). We attribute a significant increase in cloud reflectivity to enhanced cloud droplet number concentrations within a major shipping corridor in the southeast Atlantic. Prevailing winds constrain emissions around the corridor, which cuts through a climatically important region of expansive low cloud cover. We use universal kriging, a classic geostatistical method, to estimate what cloud properties would have been in the absence of shipping. In the morning, cloud brightening is consistent with changes in microphysics alone, whereas in the afternoon, increases in cloud brightness from microphysical changes are offset by decreases in the total amount of cloud water. We calculate an effective radiative forcing within the southeast Atlantic shipping corridor of approximately −2 W/m2. Several years of data are required to identify a clear signal.


Tracks Without a Trace?
For decades, ships burning high-sulfur-content fuels have crossed the world's oceans, emitting airborne particles (aerosol) and aerosol-precursor gases in regions with relatively low levels of natural aerosol (Capaldo et al., 1999;Eyring et al., 2010). Changing cloud reflectivity due to interactions with aerosol particles has long been a major driver of uncertainty in assessments of present and future anthropogenic impacts on Earth's climate (Andreae et al., 2005;Myhre et al., 2013). Marine clouds, such as the subtropical stratocumulus decks that cover much of the low-latitude oceans, are particularly sensitive to increases in aerosol concentrations as background aerosol levels are much lower over the oceans than over land (Oreopoulos & Platnick, 2008). To date, however, studies using satellite remote sensing have been unable to determine whether shipping emissions have a discernable, climatically relevant impact on cloud-radiative properties (Peters et al., 2011;Schreier et al., 2007). This is surprising, as global climate models indicate cloud responses to shipping emissions should produce a substantial radiative effect (Capaldo et al., 1999;Lauer et al., 2007;Partanen et al., 2013;Peters et al., 2014;Righi et al., 2011;Sofiev et al., 2018).
Ship tracks, or trails of cloud perturbations associated with emissions from individual ships, have been studied since the mid-1960s. Multiple hypotheses, such as that the tracks were aircraft contrails or even secret missile tests, were considered before they were correctly identified as resulting from ships traveling through conditions of shallow, cloudy marine boundary layers (MBLs) with low background aerosol levels (Conover, 1966;Twomey et al., 1968). As aerosol concentrations increase-as happens when ships emit carbonaceous particles directly and the sulfur dioxide (SO 2 ) produced by burning shipping fuel is oxidized to create sulfate (SO 4 )-more cloud condensation nuclei (CCN) are available to form liquid cloud droplets (Capaldo et al., 1999;Hobbs et al., 2000). Assuming the amount of liquid water in the clouds remains constant, this increases the cloud droplet number concentration (N d ) and decreases the effective radius (r e ) of the droplets, resulting in more reflective clouds in what is known as the effect (Twomey, 1974(Twomey, , 1977. In the 1980s and 1990s, satellite (Coakley et al., 1987) and aircraft (Durkee, Chartier, et al., 2000;Durkee, Noone, et al., 2000;Radke et al., 1989) measurements confirmed that cloud condensation nuclei from shipping emissions increases N d and decreases r e , leading to the expected increase in cloud optical thickness (τ), a measure of cloud brightness.
Rapid adjustments in cloud macrophysics can either enhance or counteract the microphysical effect on cloud brightness in warm (liquid phase) clouds. Initial ship track studies seemed to suggest that the liquid water path (L), or total amount of condensate, increased in ship tracks because shifting the droplet size distribution to smaller radii can reduce the loss of water from drizzle (Albrecht, 1989). However, later analyses (Ackerman et al., 2000;Chen et al., 2012;Coakley & Walsh, 2002) demonstrated that this is not always the case: Cloudiness can also be reduced, as increasing the cloud droplet number concentration enhances the cloud top entrainment of dry air, drying and deepening the MBL (Ackerman et al., 2004;Bretherton et al., 2007;Seifert et al., 2015;Wood, 2007). There is observational and model evidence of higher cloud tops in some ship tracks as compared to the surrounding cloud fields, which could be related to both of the aforementioned adjustments (Christensen & Stephens, 2011;Taylor & Ackerman, 1999).
In this paper, we focus on aerosol-cloud interactions in liquid clouds only, although it should be noted that shipping emissions may also be important for mixed-phase and ice cloud properties (Christensen et al., 2014;Possner et al., 2017). Indeed, lightning appears to be enhanced over major shipping corridors in the northeastern Indian Ocean and the South China Sea, which has been hypothesized to be due to convective invigoration from shipping-related aerosol perturbations (Thornton et al., 2017).
Although ship tracks have provided invaluable testbeds for aerosol-cloud interaction (ACI) hypotheses, observations have thus far suggested their overall climatic importance to be rather limited (Schreier et al., 2007). This may be due to undercounting, as clearly visible ship tracks are relatively rare (Gryspeerdt, Smith, et al., 2019). Hundreds or thousands (Campmany et al., 2009;Christensen et al., 2014;Gryspeerdt, Goren et al., 2019;Toll et al., 2017Toll et al., , 2019 have been identified by satellite per year even though there are on order 100,000 ships in the global fleet (Eyring et al., 2010). However, high-resolution model results suggest there can be radiatively important impacts even when a clear ship track is not easily discernible (Possner et al., 2018). In an attempt to capture shipping's effect more holistically, one study looked for noticeable changes in cloud properties upwind or downwind of shipping corridors by following near-surface air mass trajectories that crossed these corridors, yet no unambiguous cloud microphysical or macrophysical changes were detected (Peters et al., 2011).
Whereas observations show a small or unclear impact of shipping on global ACI, climate models (Capaldo et al., 1999;Lauer et al., 2007;Partanen et al., 2013;Peters et al., 2014;Righi et al., 2011;Sofiev et al., 2018) have produced substantial radiative forcings (energy imbalances) ranging from −0.06 to −0.6 W/m 2 . A model known to produce sizable ACI effects from shipping emissions was sampled using the aforementioned trajectory technique of Peters et al. (2011) but also could not produce a clear signal, demonstrating that detection is hampered by the high meteorological variability of tropical low cloud properties and the lack of knowledge of what cloud properties would be in the absence of anthropogenic aerosol perturbations (Peters et al., 2014).
If the shipping signal is so elusive at regional and global scales, why not simply focus our efforts elsewhere? Numerous observational studies have found a correspondence between aerosol concentrations and cloud properties (Bréon et al., 2002;Martin et al., 1994;Nakajima et al., 2001;Quaas et al., 2008;Sekiguchi et al., 2003). However, the confounding of aerosol effects and other meteorological variations that can significantly influence cloud properties remains a serious challenge for disentangling the magnitude of the aerosol effects alone (Adebiyi et al., 2015;Gryspeerdt et al., 2016;Stevens & Feingold, 2009). For this reason, "natural experiments" in which there is a clear aerosol perturbation independent of meteorological influence, such as volcanic eruptions (Gryspeerdt, Goren et al., 2019;Malavelle et al., 2017;McCoy & Hartmann, 2015;Toll et al., 2017Toll et al., , 2019 and ship tracks (Chen et al., 2012;Gryspeerdt, Goren et al., 2019;Toll et al., 2017Toll et al., , 2019, may represent our best opportunity for constraining ACI absent controlled experiments .

General Approach
A major shipping corridor in the southeast Atlantic was excluded from the previously discussed analyses of Peters et al. (2011Peters et al. ( , 2014 because the vast majority of their MBL air mass trajectories run parallel to the shipping corridor as opposed to crossing it. Indeed, MBL winds blow almost perfectly parallel to the shipping corridor as it cuts through an extensive stratocumulus cloud deck (Figure 1). These winds keep the emissions confined to a relatively narrow region surrounding the shipping corridor. Here, we take advantage of this confinement to estimate "counterfactual" fields (i.e., what would be expected in the absence of shipping perturbations) for the shipping corridor. The counterfactual fields are estimated with data from nearby, non-shipping-affected locations and covariate information using universal kriging (Zimmerman & Stein, 2010), a classical geostatistical method for spatial interpolation (see section 2.4 below). The counterfactual fields can then be compared to the "factual" fields of reanalysis or satellite data that are believed to be affected by shipping. This is conceptually akin to running a climate model with emissions turned off or on, and we refer to the kriged counterfactual fields as "NoShip" and the observation/reanalysis factual fields as "Ship" in an extension of this analogy. It should be emphasized that this method analyzes the shipping corridor "top-down" as a whole using climatological fields rather than "bottom-up" via the aggregation of individual ship tracks.
We restrict our analysis to the 2003-2015 climatology (as that is the maximum time frame for which all of our data sources are available) and to the austral spring (September-October-November) season in which both cloud coverage and shipping emissions maximize in the region. Because the shipping lane spans a Figure 1. Austral spring shipping emissions, meteorology, and cloud properties in the southeast Atlantic. (a) SO 2 emissions flux from international shipping (shading) and reanalysis winds at 1,000 hPa (barbs; half line = 5 m/s, full line = 10 m/s, winds blow from tail to head) for austral spring (September-October-November). (b) Satellite-derived cloud droplet number concentration (shading) and cloud fraction (contours of 80% and 90%) for austral spring. White boxes mark the tropical and subtropical regions of analysis.

Data
All variables analyzed in this study using the approach described in section 2.1 are listed in Tables 1  and 2. We use cloud property retrievals from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument onboard the Terra and Aqua satellites, radiative fluxes from the Clouds and the Earth's Radiant Energy System (CERES) multi-sensor satellite products, and meteorological and aerosol properties from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2). SO 2 emissions from international shipping are taken from the Emissions Database for Global Atmospheric Research (EDGAR) for September, October, and November 2010 (Crippa et al., 2018). Sector-specific data broken down by month is only available for 2010. Because shipping emissions in the southeast Atlantic do not appear to undergo any major trends (Crippa et al., 2018), the 2010 values should be representative of all years in this study. SO 2 data are interpolated from the native 0.1°× 0.1°resolution to the 1.0°× 1.0°g rid used for the monthly average MODIS products discussed below.
Meteorological (1,000 hPa horizontal winds and potential temperatures at 800 and 1,000 hPa) and aerosol (surface sulfate and black carbon mass concentrations) data in this study come from MERRA-2 (Gelaro et al., 2017;Randles et al., 2017). Shipping emissions from the EDGAR database are incorporated into the Goddard Chemistry Aerosol Radiation and Transport model, which serves as the aerosol module for MERRA-2. MERRA-2 data are interpolated from the native 0.5°× 0.625°resolution to the 1.0°× 1.0°M ODIS grid.
Cloud fraction ( C ), cloud effective radius (liquid phase), cloud optical thickness (liquid phase), and liquid water path are taken from the monthly Level-3, Collection-6 MODIS instrument products (Hubanks et al., 2019) for both Terra (daytime satellite overpass time~10:30 local and nighttime satellite overpass time~22:30 local) and Aqua (daytime satellite overpass time~13:30 local and nighttime satellite overpass time~01:30 local). Except for the cloud fraction fields, MODIS cloud properties are only available for the daytime overpass times. Cloud droplet number concentration is calculated using the retrievals of effective radius, cloud optical thickness, and cloud top temperature from MODIS/Aqua assuming a subadiabatic "Idealized Stratiform Boundary Layer Cloud" model (Bennartz & Rausch, 2017).
Daily average values of top-of-atmosphere (TOA) all-sky and clear-sky shortwave and net radiative fluxes and cloud area fraction are taken from the CERES Energy Balanced and Filled (EBAF) top-of-atmosphere Edition-4.0 data product (Kato et al., 2018;Loeb et al., 2018 Two additional data sources are utilized to calculate observationally informed estimates of the global radiative forcing (RF) due to the Twomey effect alone (R Twomey ), cloud adjustments alone (R adj ), and the total effective radiative forcing (ERF ACI , or E when used in equations) due to aerosol-cloud interactions including cloud adjustments in low clouds. Surface sulfate mass concentration data are analyzed from "historical" runs of those global climate models participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6) that had posted relevant data as of January 20, 2020. More information about the models used and their ensemble members is provided in supporting information Table S1. The average sulfate mass concentration of the lowest model level for the first 15 years of the historical run (January 1850 to December 1864) is taken as the "pre-industrial" value, and the corresponding average for the last 15 years (January 2000 to December 2014) is taken as the "present-day" value. Although sulfur dioxide from certain industries may be injected above the lowest model level in some models, it should be mixed efficiently in the boundary layer, so the monthly sulfate mass values used should be representative of sulfate throughout the boundary layer.
Radiative forcing estimates are calculated separately for "cumuliform-type" clouds (for which we consider our results from the tropical region to be representative) and "stratiform-type" clouds (for which we consider our results from the subtropical region to be representative). The Extended Edited Cloud Reports Archive (EECRA), which uses visual observer data from weather stations and ships that is reported in the World Meteorological Organization's (WMO) synoptic code (Eastman et al., 2011;Warren et al., 1986Warren et al., , 1988, is used to create cumuliform and stratiform cloud fractions. The cumuliform cloud fraction is created by aggregating observations of Types 1-2 (cumulus) and Type 8 (cumulus under stratocumulus) whereas the stratiform cloud fraction is created by aggregating observations of Types 4-5 (stratocumulus), Type 6 (stratus), and Type 11 (fog) from between 06:00 and 18:00 local time for 1954-2008 (ocean observations) and 1971-2009 (land observations). Type 11 is not contained in the original WMO synoptic code and is processed specifically for the EECRA data set. The exclusion of Type 3 (cumulonimbus without anvils), Type 7 (stratus and cumulus fractus of bad weather), and Type 9 Note. All values reported are averaged over the grid boxes with a relative increase in SO 4 exceeding 20%. 95% confidence intervals are indicated in parentheses. Variables that are both field significant at the 95% confidence level and have average Ship-NoShip differences that are significantly distinguishable from zero at the 95% confidence level are in bold.
(cumulonimbus with anvils) from our analysis should eliminate almost all observations coinciding with heavy precipitation.
The CMIP6 and EECRA data are interpolated or aggregated, respectively, to a common 5.0°× 5.0°grid. We require 20 valid observations per grid box to compute an average for the EECRA cloud fractions. Land boxes with multiple stations are averaged by first computing mean statistics for each station within the box; those means are then averaged (no weighting applied). Averages for boxes containing both land and ocean are calculated by computing an average for the land and ocean portions separately. Those separate averages are then weighted by land/ocean fraction in each grid box. Missing data are filled by setting the grid box average to zero for areas poleward of 67°(this mainly affects Greenland and Antarctica) and by averaging all neighboring valid grid box values otherwise (this mainly affects certain subtropical desert regions like the Sahara).

Albedo Calculations
Cloudy-sky, or overcast, albedo (A cld ) is calculated from the total scene albedo (A), clear-sky scene albedo (A clr ), and cloud fraction (C) by rearranging the equation: Note. All values reported are averaged over the grid boxes with a relative increase in SO 4 exceeding 20%. 95% confidence intervals are indicated in parentheses. Variables that are both field significant at the 95% confidence level and have average Ship-NoShip differences that are significantly distinguishable from zero at the 95% confidence level are in bold.
The overcast albedo calculated using equation (1) represents the total scene albedo when clouds are present, not the albedo of the clouds themselves, as the overlying atmosphere also absorbs and scatters light.
(Hereafter, scene albedos, or the albedos as seen from space, are denoted as A whereas the albedos of individual constituents are denoted as α.) Following a simplified single-layer atmospheric model for solar radiation (Donohoe & Battisti, 2011;Qu & Hall, 2005), we estimate the albedo of the atmosphere (α atm ) when clouds are not present by numerically solving the equation: where α sfc is the surface albedo and T is the transmissivity of the atmosphere (calculated as the ratio of downwelling shortwave radiation reaching the surface to that incident at TOA). The denominator of the second term on the right-hand side of the equation accounts for multiple scattering off the atmosphere and surface.
To calculate cloud albedo, we modify equation (2) under the assumption that in overcast scenes with MBL clouds, reflection from the surface makes a negligible contribution to the total outgoing shortwave flux.
Instead of considering the full atmospheric column, we now want to isolate the extinction of light in the free troposphere (FT) above the cloudy boundary layer. We estimate the albedo of the free troposphere (α FT ) by scaling the atmospheric albedo from equation (2) by the ratio of extinction occurring in the free troposphere to that for the full atmospheric column: where T FT is the transmissivity of the free troposphere (technically, the full atmospheric column minus the boundary layer), as calculated using the CERES SYN downward shortwave fluxes at the TOA and 850 hPa levels.
Cloud albedo (α cld ) is then calculated by numerically solving the equation: Conversions can be made between changes in overcast albedo and changes in cloud albedo by differentiating equation (4): In regions like the southeast Atlantic where few high clouds are present, the interpretation of α cld as calculated via equation (5) is more straightforward than in other locations where cloud types with differing characteristic albedos vary greatly, such as the midlatitudes. For global values of cloud albedo in liquid phase clouds, we instead calculate α cld using MODIS liquid cloud optical thickness under the Eddington approximation for non-absorbing media (e.g., Segrin et al., 2007): where the asymmetry parameter g is assumed to be 0.85.

Universal Kriging
Universal kriging is a classic geostatistical method (Zimmerman & Stein, 2010) designed to estimate some value at unknown spatial locations based partially on nearby observations of the same value. At each unknown location, the mean value is estimated from a regression model. Error, or noise, around all mean values is assumed to be spatially correlated. The correlation of the error between two values is further assumed to be a function only of the distance between locations, a property known as stationarity. We use the statistical package geoR (Ribeiro & Diggle, 2018) in R (R Core Team, 2019) for implementing this analysis.
Shipping-influenced grid boxes are identified by choosing the 1.0°× 1.0°grid box at each latitude (within the longitude ranges specified above for each region) with the maximum SO 2 emissions from the EDGAR emissions database and the two neighboring grid boxes to the east and west. This results in five grid boxes identified as "shipping-affected" for each latitude in the two analysis boxes, for a total of 40 grid boxes per region. This extent was chosen to ensure that most of the area affected by the diffuse edges of the aerosol perturbation were categorized as shipping-influenced, while maintaining enough reference grid boxes to robustly fit the kriging model.
The mean value is obtained with linear regression using the possible covariates of latitude (lat), longitude (lon), their squares (lat 2 and lon 2 ), and their product (lat*lon) in addition to lower tropospheric stability (Δθ; defined here as the potential temperature difference between 800 and 1,000 hPa) and an "effective" measure of LTS, (Δθ) eff , accounting for horizontal MBL advection: where the zonal and meridional winds (u ! ) are taken from the 1,000 hPa level. This second measure was created to account for the fact that cloud cover is more strongly correlated with the lower tropospheric stability the MBL experienced 24-36 hr prior than with its instantaneous value (Klein et al., 1995;Mauger & Norris, 2010)-this effective lower tropospheric stability measure is essentially a Lagrangian adjustment to what is otherwise a fundamentally Eulerian analysis. Not all potential regressors are used to create each variable's mean function. To select an appropriate model, the Bayesian information criterion (BIC) is computed for all possible combinations of regressors and the combination that minimizes the BIC is selected. Tables S2 and S3 list the regressors used for each variable for the subtropical and tropical regions, respectively.
To ensure errors around the mean function are normally distributed (insofar as possible), some variables are transformed via the logarithm or logit function. Each variable's transformation, if any, is reported in Tables S2 and S3.
The stationary error term is estimated by fitting a parametric covariance model (here assumed to be exponential) to an empirical variogram, a plot of the squared difference between pairs of variables versus their distance, using weighted least squares. Binned empirical variograms and the fitted variograms are provided for the subtropical and tropical domains in Figures S1 and S2, respectively.

Significance Testing
We perform three distinct tests of significance in this study: (1) the significance of individual grid boxes, (2) the significance of the field of grid boxes, and (3) the significance of the average Ship-NoShip difference within the core shipping lane.
For the first test, each grid box is considered to be individually significant if its factual (Ship) value is either above the 97.5th or below the 2.5th percentile of the distribution obtained via kriging for its counterfactual (NoShip) value. For the second test, we evaluate how extreme the number of individually significant grid boxes in the full region is compared to the number we would expect under the null hypothesis that the region is unaffected by shipping. Using the statistical model that kriging provides, we simulate 5,000 null fields for the full region. The p value is the fraction of simulations that have a number of individually significant grid boxes equal to or greater than that of the factual case. By simulating full regions in addition to individual grid boxes, we account for the effect correlated error structures can have. Because we are testing multiple hypotheses, we apply a Benjamini-Hochberg adjustment (Benjamini & Hochberg, 1995;Ventura et al., 2004) to this p value to control the false discovery rate, or the proportion of false positives expected given the number of tests. Accounting for multiple testing meaningfully affects our results: Both the CERES EBAF total albedo and MODIS/Aqua liquid water path Ship-NoShip differences in the subtropical region would be considered field significant (at the standard 5% significance level) using the raw p values but are not significant under the more stringent adjustment.
If a variable shows individual grid boxes both above and below the 95% confidence interval, field significance is assessed using the direction in which the majority of the individually significant grid boxes fall. The only variable that is assessed by our method to be field significant in which this situation arises is the MODIS/Aqua 13:30 local cloud fraction in the tropical domain, in which seven grid boxes lie above the 97.5th percentile and two grid boxes lie below the 2.5th percentile. Although the majority of the individually significant grid boxes show an increase in cloud fraction, the central tendency of the field is toward reductions in cloud fraction. While it is possible that there are real cloud responses in this region that are subtler than our methods were designed to accommodate, we consider it more likely that this is a false positive.
For testing if the average Ship-NoShip difference is significant, we focus on the core of the shipping corridor, defined as grid boxes with MERRA-2 SO 4 Ship-NoShip increases greater than 20%, rather than the more diffuse edges of the aerosol anomaly. We deem the average change to be significant if the resulting values are distinct from zero at the 95% confidence level. We estimate uncertainty as the spread in the factual-counterfactual (Ship-NoShip) differences for all 5,000 simulations from the statistical model obtained from kriging.
We place the greatest amount of confidence in the results that are both field significant and with average Ship-NoShip difference distinguishable from zero, such as those for cloud droplet number concentration, effective radius, and cloud albedo in the subtropical region. Although some results-like those for cloud optical thickness from MODIS/Terra and liquid water path from MODIS/Aqua in the subtropical region -are not field significant, we still have moderate confidence in their reliability given the consistency of the sign of their effect (i.e., the Ship-NoShip differences within the core corridor are statistically distinguishable from zero). We have low confidence that variables that are neither field significant nor with average Ship-NoShip difference distinguishable from zero are perturbed by the shipping effects and thus refrain from drawing strong conclusions about changes in these variables.

Cloud Microphysical Changes
Tables 1 and 2 contain our results for each variable analyzed for the subtropical and tropical regions, respectively. Rows containing variables in which the Ship-NoShip values are both field significant and distinguishable from zero at the 95% confidence level have been bolded. In both regions, MERRA-2 surface black carbon (BC) and sulfate (SO 4 ) mass concentrations are statistically significant on both counts. EDGAR emissions are an input for MERRA-2, so this is not entirely surprising, although it does establish that our method is sensitive enough to distinguish a shipping perturbation above the background aerosol concentrations. Figure 2 shows the mean counterfactual NoShip estimate obtained from kriging and the factual Ship fields for three aerosol and cloud microphysical variables: surface sulfate mass concentrations, cloud droplet number concentration, and effective radius (from MODIS/Aqua). Substantial differences between the observed Ship and mean NoShip estimate indicate an effect of shipping on the cloud properties.
Although the surface sulfate mass increases significantly in both the tropical and subtropical regions, only the latter region shows consistently significant changes in cloud microphysical properties. From Tables 1 and 2, we can further see that only two satellite-based variables (MODIS/Aqua N d and CERES SYN total albedo in the morning) pass both significance tests in the tropical region (and even then, with p values well above 0.01), whereas nine (MODIS/Terra r e , MODIS/Aqua N d and r e , CERES EBAF overcast albedo and net flux, CERES SYN-M cloud albedo, and CERES SYN-MHour total albedo, overcast albedo, and cloud albedo in the morning) do in the subtropical region (with p values all well below 0.01). We therefore focus the bulk of our discussion on the subtropical region.
A commonly used metric (Feingold et al., 2001;McComiskey et al., 2009;McCoy et al., 2017) for the strength of aerosol-cloud interactions relates the relative change in cloud properties with the relative change in CCN concentrations (or some suitable proxy, like aerosol optical depth or sulfate mass). Here, we define such an "ACI parameter" (also commonly known as an "Indirect Effect parameter"), β, as the fractional change in N d per fractional change in sulfate mass: In the subtropical region, an approximately 25% increase in SO 4 corresponds to an approximately 5% increase in N d , resulting in an ACI parameter of 0.21 (95% confidence interval: 0.16 to 0.27).
Estimating the ACI parameter with r e instead of N d , δln Nd ½ ð Þ , results in higher values of 0.31 (0.22 to 0.40) and 0.32 (0.25 to 0.40) for MODIS/Aqua and MODIS/Terra, respectively. However, this calculation assumes that liquid water path does not change, which is not necessarily true. The discrepancy between ACI parameters estimated using N d versus r e is a cause for caution in interpreting r e -derived ACI values when the liquid water path may also be changing. It should also be noted that the susceptibilities we derive using climatological fields are valid for O(100 km) spatial scales and decadal timescales. Care should be taken when comparing to values derived at finer resolution, for example, via instantaneous aircraft measurements.
Although no published Terra N d product is currently available, the relative change in N d can be estimated from the changes in r e and liquid water path as From equation (9), δN d /N d is estimated as 7.2% (5.2 to 9.3%) in the morning. The ACI parameter calculated using the MODIS/Terra δN d /N d is 0.30 (0.21 to 0.39). One physical explanation for larger cloud droplet number increases in the morning is that higher in cloud supersaturations from stronger updrafts in the morning may activate a larger fraction of the shipping-generated aerosol particles, which tend to be quite small in diameter Petzold et al., 2010). The MBL is also more likely to be decoupled in the afternoon, which could inhibit transport of aerosol from the ships to cloud base.

Cloud Macrophysical Changes
Despite the clear increases in N d , we do not find significant and consistent changes in cloud fraction in the daily average or at any particular time of day from either MODIS or CERES satellite retrieval products (Tables 1 and 2). We therefore refrain from drawing firm conclusions about cloud fraction changes due to shipping in either region besides that large changes in either direction are unlikely. More modest increases or decreases cannot be ruled out.
In the subtropical region, we find a decrease in liquid water path in the afternoon that is clearly distinguishable from zero within the estimated uncertainty. The central tendency of liquid water path is a mild decrease during the morning, although no change or a small increase are also plausible within the uncertainty. Liquid water path changes in the tropical region are highly uncertain in both magnitude and sign.
The larger decrease in liquid water path as the day progresses in the subtropical region is consistent with previous modeling results (Sandu et al., 2008) showing that the net sign of the cloud adjustments at night is sensitive to the competition between drizzle suppression and entrainment drying whereas during the day the enhanced entrainment effects dominate. This result conflicts, however, with the findings of other observational studies from manually detected ship tracks in the northeast Pacific that show a larger liquid water path decrease for MODIS/Terra than MODIS/Aqua (Christensen et al., 2009;Segrin et al., 2007). It is unclear whether this reflects a true difference in cloud adjustments between regions or rather reflects differences in the statistics of individual ship tracks versus the aggregated signal (including both clearly detectable and subtler shipping influences) within the shipping corridor.
We can define a susceptibility parameter, λ, relating the relative change in liquid water path to a relative change in cloud droplet number concentration: For MODIS/Terra, λ is −0.09 (95% confidence interval: −0.31 to 0.08) and for MODIS/Aqua, λ is −0.39 (−0.68 to −0.14), giving a daily mean λ value of −0.24 (−0.42 to −0.08). These λ values are at the more negative end of estimates previously reported in the literature, as compiled by Sato and Suzuki (2019), and are in line with the values for higher N d from Gryspeerdt, Goren et al. (2019) and lower r e from Toll et al. (2019). However, as the relationship between liquid water path and cloud droplet number is not necessarily linear in relative changes or even monotonic (Gryspeerdt, Goren et al., 2019), there is uncertainty inherent in applying the λ values estimated here more broadly (particularly in regions with greater precipitation).
Heating of the marine boundary layer (a "semi-direct" aerosol radiative effect) by black carbon, which is also emitted by ships, could also potentially lead to decreases in liquid water path during the day. A diabatic heating rate is calculated for the absorption of black carbon throughout the depth of the marine boundary layer in the afternoon for the subtropical domain. An absorption coefficient (k a ) for the black carbon emitted by shipping is calculated by assuming a relatively high value for the mass absorption cross section of black carbon of 15 m 2 /g (Bond et al., 2013) and multiplying by the estimated increase in black carbon from Table 1 and the air density (ρ a ). The amount of incoming energy flux absorbed (dF) in an infinitesimal vertical layer (dz) of boundary layer air is given by Beer's Law as Integrating equation (11) over the depth of the marine boundary layer (h) and recognizing that k a h ≪ 1, the shortwave flux absorbed in the marine boundary layer (ΔF) can be expressed as where F ⊙ is the top-of-atmosphere solar shortwave flux, here averaged around the Aqua overpass time over the subtropical domain (~1,080 W/m 2 ). Using the top-of-atmosphere insolation overestimates the shortwave flux available to be absorbed by neglecting scattering and absorption of sunlight by gases and particles above the MBL. The diabatic heating rate is then calculated as where T is temperature, t is time, and c p is the specific heat capacity of dry air. For our estimated values, the heating rate in the afternoon from equation (13) is only~0.001 K/hr, which is negligible. Semi-direct effects, therefore, are unlikely to be important for the cloud adjustments we observe.
Cloud optical thickness is related to the effective droplet radius and liquid water path as Figure 3 decomposes the relative cloud optical thickness changes into components due to relative changes in liquid water path and r e separately for the morning (Terra) and afternoon (Aqua): Probability densities are calculated from the 5,000 Ship-NoShip estimates for each variable via Gaussian kernel density estimation (KDE). The relative cloud optical thickness change calculated from the liquid water path and r e components closely matches the relative change calculated directly from the Ship-NoShip differences.
Although the increase in cloud optical thickness is dominated by the decrease in r e in the morning, by afternoon, the decrease in liquid water path is sufficiently large to offset this effect. At least in this region, cloud adjustments partially counteract rather than enhance the Twomey effect.

Radiative Impact
The estimated change in daily average total scene albedo for the subtropical region, given the mean austral spring incoming solar radiation, leads to an effective radiative forcing of −1.9 W/m 2 (95% confidence interval: −2.9 to −0.8 W/m 2 ) in the shipping corridor (Figure 4). (The negative sign of the radiative forcing indicates that more shortwave energy is now leaving the Earth system.) Estimating the radiative forcing via the CERES EBAF net TOA flux directly leads to a somewhat more negative estimate of −2.2 W/m 2 (−3.2 to −1.2 W/m 2 ). The absolute values of these radiative forcing estimates are more than an order of magnitude greater than previous estimates for the southeast Atlantic region obtained from clearly visible ship tracks alone (Schreier et al., 2007).
The total albedo changes can be broken down into components due to changing cloud fraction and those due to changing cloud brightness as We find that changes in cloud albedo alone explain virtually all of the observed change in total albedo (Figure 4). The radiative forcing from cloud fraction changes alone is small compared to that from changing cloud albedo and is not distinguishable from zero at the 95% confidence level. Hence, it follows that changes in cloud brightness, not fractional coverage, drive the radiative response to shipping in this region.
Using the MODIS and CERES SYN-MHour data, we can analyze how the radiative forcing varies diurnally and break down the effective radiative forcing into components related to the Twomey effect and cloud adjustments. We estimate the radiative forcing due to the Twomey effect alone (Platnick & Twomey, 1994;Quaas et al., 2008) as and the radiative forcing due to cloud adjustments alone (Quaas et al., 2008) as where we use equation (5) to convert between cloud albedo and overcast (scene) albedo and have neglected adjustments due to cloud amount changes, as we cannot estimate the sign and magnitude of potential cloud fraction changes inside the shipping corridor with any confidence. We estimate the effective radiative forcing from increasing cloud brightness from the MODIS (E MODIS ) and CERES (E CERES ) data separately as Figure 5 shows the radiative forcing estimates calculated for the Terra and Aqua overpass times. Note that for the Terra overpass time, the relative change in N d is inferred from equation (9). Interestingly, the effective radiative forcing from the CERES data is systematically larger than that calculated from MODIS, although they are still consistent within the estimated uncertainty. In the morning, cloud brightening from the Twomey effect is consistent with the main area of overlap between the CERES and MODIS effective radiative forcing estimates, whereas during the afternoon, the radiative forcing from the Twomey effect is larger than the main area of overlap between the CERES and MODIS effective radiative forcing estimates. This is consistent with the decrease in liquid water path in the afternoon counteracting much of the brightening from increasing the cloud droplet number concentration. In the morning, when complications due to changes in liquid water path are expected to be small, the general agreement between the radiative forcing from the Twomey effect and the total effective radiative forcing estimates from cloud brightening (particularly E MODIS ) suggests that the Platnick and Twomey (1994) method of turning a relative N d change into a radiative forcing works reasonably well.
The central estimates of the cloud albedo increases in the tropical region are also broadly consistent with a dominant Twomey effect. However, because the cloud albedo changes are neither field significant nor statistically distinguishable from zero at the 95% confidence level, we refrain from drawing strong conclusions about cloud brightening in the tropical region.

Timescales of Detectability
Here, we were able to demonstrate for the first time that detectable changes in cloud radiative forcing are generated by commercial shipping on climate-relevant regional spatial scales and decadal temporal scales. However, we were only able to detect a robust response in the stratocumulus-dominated subtropical domain, not in the tropical domain where cumuliform clouds are more common. We argue that this may be a signal-to-noise question. In regions where the sensitivity of cloud properties and organization to meteorological variability is larger, such as the more tropical domain (George & Wood, 2010;Peters et al., 2014), longer time series may be needed to average out the variability and to obtain the signal. Previous attempts to detect shipping perturbations were based on merely a few years of data, which may not have been sufficient to smooth out the noise. With only 1 or 2 years of data, it is possible to surmise the sign of the N d effect but not its magnitude ( Figure 6). Neither are necessarily clear for the A cld effect. The estimates better approximate their climatological distributions once 5 to 6 years of data are available.

To test this hypothesis, we reprocess our results for MODIS/Aqua N d and CERES EBAF
Given the amount of time needed to detect clear signals for the relatively uniform stratocumulus clouds in the subtropical domain, the hypothesis that tropical clouds require a substantially longer record to be able to detect a signal above the noise seems plausible. As a further complication, remotely sensed cloud properties are more uncertain in the tropical region due to the lower frequency of occurrence of clouds, three-dimensional radiative effects, heterogeneity of the cloud field, cloud-edge effects, and other biases (Bennartz & Rausch, 2017;Grosvenor et al., 2018).
It is also, of course, possible that the lack of detection of a strong shipping signal is not a result of observational limitations but rather a reflection that shipping effects truly are smaller in the tropical region. For instance, the marine boundary layer is deeper and more likely to be decoupled in this region, so the aerosols and precursor gases emitted by ships near the surface may not be efficiently transported to the subcloud layer. Additionally, it is possible that cumuliform clouds are less susceptible to aerosol-driven increased in N d than stratiform clouds for physical reasons, such as differences in supersaturation adjustments between the two regimes (Jia et al., 2019).

Implications for 2020 International Maritime Organization Regulations
Our findings are undergoing a real-world test as International Maritime Organization (IMO) regulations (International Maritime Organization, 2016;Sofiev et al., 2018) limiting fuel sulfur content to 0.5% by mass (from the present standard of 3.5%) came into effect on 1 January 2020. Although the fuel used will still contain some sulfur content, the particles produced by ships under the new regulations may be too small to effectively serve as cloud condensation nuclei at the relatively low supersaturations characteristic of marine stratocumulus clouds (Petzold et al., 2010). An analysis of the effect of reducing fuel-sulfur content restrictions from 1.0 to 0.1% in 2015 within the emission control area off the western coast of North America found a dramatic decrease in manually detected ship tracks immediately following the strengthening of the regulation (Gryspeerdt, Smith, et al., 2019). The 2020 global IMO regulations are therefore likely to dramatically reduce the cooling effect of aerosol-cloud interactions associated with the global shipping industry. If our analysis is valid, we should be able to observe decreases in N d and overcast albedo by the mid-2020s assuming the regulation achieves broad compliance and is adequately enforced.

Implications for Marine Cloud Brightening
There are also practical implications of our detectability findings for the study of aerosol-cloud interactions generally, and in particular for possible field tests of marine cloud brightening as a geoengineering response to anthropogenic climate change (Latham et al., 2012;Wood & Ackerman, 2013;Wood et al., 2017). The first paper documenting ship tracks in the 1960s remarked that it may be possible to cool the climate via cloud-seeding over the oceans (Conover, 1966). More recently, interest has grown in exploring the feasibility of doing so to counteract warming due to anthropogenic greenhouse gas emissions (Latham, 1990;Latham et al., 2012;Lawrence et al., 2018). It has even been proposed that the international shipping fleet could alter the fuel it burns from low-sulfur near the coasts to higher-sulfur further offshore to capture some of the health benefits of reduced air pollution while maintaining most of the ships' cooling effects (Partanen et al., 2013). The results here suggest that a regional-scale field test of marine cloud brightening in stratocumulus-dominated regions would likely be successful in terms of increasing cloud albedo; however, it may require several years for the regional perturbation to the outgoing shortwave energy flux to be clearly detectable via satellite remote sensing. This agrees well with previous work suggesting that multiple years of observations may be required to detect deliberate albedo modifications depending on the size and abruptness of the perturbation (Seidel et al., 2014). If detection from space is considered a key criterion of success for the field test, this relatively large time commitment raises the question of at what point a regional field test blurs into regional implementation of marine cloud brightening (Robock et al., 2010).

Effective Radiative Forcing
If we assume that that relationships between sulfate, cloud droplet number concentration, and cloud albedo calculated within the shipping corridor in the southeast Atlantic hold globally, we can estimate the present-day effective aerosol-cloud radiative forcing for a given increase in sulfate from pre-industrial times. The global estimate is obtained by combining estimates for cumuliform-type clouds using relationships calculated from the tropical domain of the shipping corridor and estimates for stratiform-type clouds using relationships calculated from the subtropical domain of the shipping corridor. Previous work has shown that the relationship between cloud droplet number concentration and MBL sulfate mass is similar across different stratocumulus-dominated ocean and land regions (McCoy et al., 2017). This suggests that our assumption that values derived for a given cloud type in one particular region can represent similar cloud types across the globe is reasonable.
The mean ratio of present-day to pre-industrial surface sulfate mass concentrations ([SO 4,PD ]/[SO 4,PI ]) from the CMIP6 models is shown in Figure 7a. The multi-model mean is calculated by first calculating the mean of all ensemble members of each model and then averaging all individual-model means.
We calculate ERF ACI by converting the ratio of present-day to pre-industrial sulfate at each grid box into an estimate of increased cloud albedo. To do this, we first define a relative cloud albedo susceptibility, χ, as and use the tropical and subtropical CERES EBAF cloud albedo and averaged MODIS/Terra and MODIS/Aqua cloud droplet number concentration changes to calculate a value of 0.22 (95% confidence interval: 0.11 to 0.33) for stratiform clouds and 0.32 (−0.13 to 1.03) for cumuliform clouds within the shipping corridor. The large uncertainty in the cumuliform cloud estimate is due to both the poor constraint on cloud albedo changes and the small magnitude of the cloud droplet number concentration changes in the tropical region. As before, the relative MODIS/Terra N d increase is calculated using equation (9).
The ratio of present-day to pre-industrial cloud albedo can then be calculated by first raising the sulfate ratio to the power of the ACI parameter to calculate the cloud droplet number concentration ratio and then raising the cloud droplet number ratio by the relative cloud albedo susceptibility: The ACI parameter ( Note that the change in sulfate is large compared to the pre-industrial sulfate concentration whereas the change in cloud albedo is generally much smaller than the pre-industrial cloud albedo. The power law relationship in equations (22) and (23), rather than a more straightforward linear multiplication of the sensitivity parameters with relative pre-industrial to present-day changes, causes the cloud droplet number concentration and cloud albedo increases to saturate more rapidly with increasing aerosol and is necessary to avoid unphysical albedo increases.
From the ratio calculated with equation (23), we estimate the absolute change in overcast albedo (ΔA cld ) as where equation (5) is used to convert between overcast and cloud albedo and the absolute value of the present-day cloud albedo (α cld ) is estimated as the average of the MODIS/Terra and MODIS/Aqua liquid cloud albedos as calculated using equation (6).
Because the EECRA stratiform and cumuliform cloud fractions are created using observations of clouds made from the surface, we must account for the masking of low cloud albedo changes by higher clouds before we can calculate the top-of-atmosphere ERF ACI . We define "high clouds" to be those with cloud top pressure values of less than 680 hPa and calculate a high cloud fraction, C high , using the joint histograms of daytime cloud fraction and cloud top pressure from the standard MODIS product (values for Aqua and Terra are averaged to create a daytime average). We assume random overlap between high and low clouds. Putting all of the above together, we calculate ERF ACI for each grid box as where the solar flux is now representative of the annual mean and C low refers to either the stratiform or cumuliform cloud fraction. ERF ACI values are calculated separately for stratiform and cumuliform clouds using the cloud fraction and albedo change associated with each cloud type and the ERF ACI for all low clouds is taken as the sum of the stratiform and cumuliform contributions. It should be noted that our estimates do not include potential cloud fraction adjustments, as we did not find clear evidence for any systematic cloud fraction changes in response to the shipping perturbation.
To generate a distribution of ERF ACI values, for each of the 5,000 Ship-NoShip estimates of β and χ, a CMIP6 model is chosen at random and then one ensemble member is chosen at random to provide the sulfate ratio that is then used to calculate the change in cloud albedo. The range of plausible ERF ACI values thus accounts for both uncertainty in the kriging model as well as uncertainty due to the inter-model spread in sulfate increases.
Figures 7b-7d provide global maps of the resulting mean ERF ACI estimates for all low clouds, stratiform low clouds only, and cumuliform low clouds only. Unsurprisingly, ERF ACI is greatest in the Northern Hemisphere, with particularly large values over Asia, North America, and the northeast Pacific Ocean. Despite having below-average increases in sulfate, the southeast Pacific and southeast Atlantic stratocumulus regions also have large ERF ACI values due to their extensive coverage of low clouds (and relative lack of high clouds).
Global average ERF ACI values are reported in Table 3 for each cloud type and their combination. The ERF ACI from stratiform clouds is approximately twice as large as that from cumuliform clouds; however, the uncertainty in the total ERF ACI from low clouds is dominated by the cumuliform cloud contribution. This is primarily due to the poor constraint from our CERES EBAF results in the tropical domain as compared to those found in the subtropical domain.
Our global mean estimate of ERF ACI in warm low clouds of −0.97 W/m 2 (95% confidence interval: −1.63 to −0.38 W/m 2 ) is constrained by causal relationships between aerosol, cloud droplet number concentration, and cloud albedo observed within the southeast Atlantic shipping corridor. As discussed in detail below, however, the forcing we calculate is broadly in line with other estimates in the literature that use regression-based analyses (in which assessing the degree of causality in observed relationships is more challenging). Our results thus represent an independent estimate that can increase our confidence in the reasonableness of the magnitude of the ERF ACI estimates based purely on regression statistics.
Our estimate is well within the ERF ACI range of −3.1 to −0.1 W/m 2 (95% confidence interval) from the global aerosol radiative forcing review resulting from the 2018 World Climate Research Programme's Grand Science Challenge on Clouds, Circulation and Climate Sensitivity workshop at Schloss Ringberg . Some discrepancy is to be expected, as potential adjustments in cloud fraction, which are not considered here, substantially lower estimates in bottom-up based global mean ERF ACI .
Assuming moderate values for the effective radiative forcing from aerosol-radiation interactions and any potential cloud fraction adjustments not accounted for in our analysis, our estimated range for ERF ACI is also consistent with global constraints placed on total aerosol radiative forcing by the historic temperature record between 1850 and 1950 as well as the hemispheric contrast in aerosol forcing. The idea to constrain global ERF ACI in this manner was first utilized in Stevens (2015) and later revised by Kretzschmar et al. (2017) and Booth et al. (2018) and places the physically plausible lower bound of global effective radiative forcing due to all aerosol effects as approximately −2 W/m 2 .
ERF ACI as assessed in the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC), −0.45 W/m 2 (90% confidence interval of −1.2 to 0.0 W/m 2 ), is less than half of our central estimate but consistent within the considerable uncertainty range (Boucher et al., 2013;Myhre et al., 2013). An analysis of climate model runs from Phase 5 of CMIP (CMIP5) found an ERF ACI estimate of −0.59 W/m 2 (2σ range of −1.21 to −0.03 W/m 2 ) when only including anthropogenic sulfate emissions and −0.92 W/m 2 (2σ range of −1.60 to −0.24 W/m 2 ) when including all aerosol sources (Zelinka et al., 2014).  As a caveat to the above, it should be noted that our estimate does not attempt to account for the different effects of aerosol-cloud interactions in mixed-phase and ice clouds, which are quite uncertain (Lohmann, 2017;Storelvmo, 2017).

Radiative Forcing Components Due to the Twomey Effect and Cloud Adjustments
We can also use our results from shipping in the southeast Atlantic to provide an observationally informed estimate of the radiative forcing from the Twomey effect alone and cloud adjustments alone through an analysis similar to the analyses found in McCoy et al. (2017) and Toll et al. (2019), respectively, and that presented above.
For the Twomey effect estimate, relative cloud droplet number changes between the present-day and pre-industrial era, Δln(N d ), are estimated using the ACI parameters (β) calculated in section 3 (defined using the relative N d increase) from MODIS/Terra and MODIS/Aqua for each region as From equation (26) From this relative droplet number change, we estimate the cloud albedo increase as where α cld is calculated from MODIS/Terra or MODIS/Aqua cloud optical depth using equation (6). This ΔA cld,Twomey value is calculated for the Terra and Aqua overpass times separately and then averaged together to create a daily mean estimate for each region/cloud type.
For the radiative forcing due to cloud adjustments (only including liquid water path changes), we calculate the relative change in liquid water path, Δln L ð Þ, by first converting the cloud droplet number ratio to a liquid water path ratio using liquid water path susceptibilities (λ) as calculated by equation (10): which we then use to calculate the change in cloud albedo due to liquid water path adjustments as As with ΔAcld,Twomey, this ΔAcld,adj value is calculated for the Terra and Aqua overpass times separately and then averaged together to create a daily mean estimate for each region/cloud type.\nThe radiative forcings due to the Twomey effect, R Twomey , and cloud adjustments, R adj , are finally calculated as where the conversion factor between cloud and overcast albedo changes is calculated for the daily mean, as in equation (24).
As with ERF ACI , results for the globally averaged R Twomey values for each cloud type and their combination are reported in  Toll et al. (2019). Our overall uncertainty for the R Twomey estimates is lower than for ERF ACI in large part because of the better constraints on MODIS cloud effective radius and cloud droplet number concentration changes as compared to changes in CERES cloud albedo in the tropical region.
Our estimate for the radiative forcing due to cloud adjustments is very small, 0.02 W/m 2 (95% confidence interval: −0.30 to 0.32 W/m 2 ), due to cancellation between competing effects from stratiform-type clouds (decreasing liquid water path) and cumuliform-type clouds (central tendency of increasing liquid water path). We have very low confidence, however, in the magnitude and even the sign of the liquid water path changes for the cumuliform-type clouds. The R adj estimate for stratiform clouds of 0.18 W/m 2 (0.05 to 0.33 W/m 2 ) is similar to the value of 0.15 W/m 2 from Toll et al. (2019) estimated using manually detected pollution tracks.
As with the radiative forcing estimates from the southeast Atlantic shipping corridor itself, the global radiative forcing estimates are larger when calculated using CERES cloud albedo as compared to MODIS-derived cloud albedo changes. It is possible that our cloud albedo estimates from section (6) are overestimated because we assume that the downwelling and upwelling transmissivity of the atmosphere is the same. In reality, absorption in certain spectral regions (like the ultraviolet ozone bands and near-infrared water vapor bands) should saturate during the first pass through the atmosphere, so the upwelling transmissivity and that relevant for any multiple scattering should be larger than that for the direct solar beam. At the same time, MODIS-derived values of cloud optical thickness and cloud droplet number concentration may be biased low due to the effects of overlying absorbing aerosols (Grosvenor et al., 2018;Meyer et al., 2015). However, our use of relative changes should blunt the effect of biases in the absolute values of the cloud properties. The cloud albedo calculated from MODIS cloud optical depth using equation (6) may also be overestimated because the Eddington approximation neglects absorption of light by cloud droplets. The absolute cloud albedo change calculated from equation (24) may therefore be biased high independent of MODIS/CERES differences. Other major sources of uncertainty in our radiative forcing estimates are addressed in section 7.4.

Spatial Patterns in Radiative Forcing
In addition to the global averages, Table 3 also records global average values for the analysis performed over the oceans and land separately and for the Northern Hemisphere and Southern Hemisphere separately. All values are reported for the globe as a whole-to get the forcing for the Northern Hemisphere alone, for example, the relevant value in Table 3 would need to be doubled.
The contribution to global radiative forcing from marine clouds is approximately twice as large as that from clouds over land. Taking into account the relative fractions of Earth's surface covered by land versus oceans, however, this implies that there has been a larger radiative forcing due to aerosol-cloud interactions over land than over the oceans since the pre-industrial era. (This is particularly true for the cumuliform cloud forcing but holds for the stratiform estimates as well.) Our analysis accounts for differences in land and ocean regions only via the differences in the coverage of stratiform and cumuliform low clouds and of high clouds, as opposed to other studies that calculate cloud property susceptibilities for land and ocean regions separately (e.g., Quaas et al., 2008;Toll et al., 2019). ACI relationships derived over land tend to show larger susceptibilities of cloud properties to changes in aerosol in cleaner as opposed to more polluted background conditions (e.g., Feingold et al., 2001), suggesting that varying background aerosol concentration may be driving the differences in susceptibility estimates between land and ocean regions. Differences in cloud susceptibility to aerosol perturbations could also be due to variation in key processes between marine and terrestrial boundary layers, however.
The Northern Hemisphere contributes three times as much forcing as the Southern Hemisphere, which is unsurprising given that the preponderance of anthropogenic pollution sources are located in the Northern Hemisphere. Table 3 further breaks down each hemisphere's forcing into ocean and land components. As would be expected, land forcing is negligible in the Southern Hemisphere, whereas in the Northern Hemisphere the land contribution is substantial. This is particularly true for the cumuliform-type cloud forcing.
We can also estimate the extent to which accounting for spatial patterns in aerosol and cloud properties matters for globally averaged estimates. A "bulk ocean" effective radiative forcing, E ð Þ ocn , is calculated by modifying equations (24)-(25) as where the ð Þ ocn overbar refers to the value being geospatially averaged over the global oceans. Depending on whether the present-day to pre-industrial sulfate ratio is calculated as the inverse of the average present-day increase or the average of the inverse of the present-day increase, E ð Þ ocn is calculated as either −1.72 W/m 2 (−2.96 to −0.73 W/m 2 ) or −0.77 W/m 2 (−1.29 to −0.34 W/m 2 ), respectively, as compared to the value of −0.62 W/m 2 (−1.08 to −0.26 W/m 2 ) from the analysis accounting for spatial covariations. Back-ofthe-envelope "bulk" calculations can thus be useful for a first-order impression of the global radiative forcing, but it is clear that the details of the spatial patterns and precise method used to generate the "bulk" quantities can matter a great deal.

Major Sources of Uncertainty
There are several sources of uncertainty that could lead to our estimates being overly strong. We implicitly assume that the lower sensitivity to aerosol perturbations generally found in continental clouds (Oreopoulos & Platnick, 2008;Quaas et al., 2008) is due to higher background aerosol concentrations in the present-day and would not necessarily apply in the pre-industrial environment. If continental regions are less sensitive to aerosol perturbations than oceanic regions more generally, applying our susceptibilities derived in a marine environment would overestimate the change over land. Similarly, the pre-industrial background state, especially for the Northern Hemisphere, is highly uncertain (Carslaw et al., 2013)-a too-pristine pre-industrial environment in the CMIP6 models would also lead to overly strong aerosol effects in our calculations. Our use of relative susceptibilities and a power law (as opposed to linear) functional form in equations (22), (23), and (28) accounts for some degree of saturation at large aerosol perturbations, although it is possible that we underestimate the degree of saturation for the very large perturbations seen over some regions like India and eastern Asia.
More broadly, assigning singular values to the cloud susceptibility parameters (β, λ, χ) implicitly assumes relationships that are monotonic and linear in relative changes. For the increase in cloud droplet number concentration with aerosol (β) and cloud albedo with droplet concentration (χ), we would expect errors here to tend to overestimate the changes (i.e., effects may saturate more rapidly in reality than under an assumption of linearity). For liquid water path changes with droplet number (λ), the effects could be subtler as factors such as the free tropospheric relative humidity and likelihood of clouds to precipitate in the absence of an aerosol perturbation could cause λ to vary in magnitude and sign with meteorological conditions (Gryspeerdt, Goren et al., 2019;Toll et al., 2017Toll et al., , 2019. The southeast Atlantic stratocumulus region during austral spring experiences very strong cloud top inversions, generally quite well-mixed conditions, and very high cloud fractions. It is quite plausible that the relationships derived in this region and season, at least for the ACI parameter (β), represent a reasonable upper-bound on what should be expected in stratiform-dominated regions elsewhere. Even our tropical, cumuliform-dominated region has quite high cloud fraction-the relatively weak relationships found in that region (which is still partially within the transition region between subtropical stratocumulus and tropical trade cumulus clouds) may be even smaller in areas that are more dominated by shallow cumuli.
Other sources of uncertainty cut in the other direction, potentially leading us to underestimate the magnitude of the radiative forcings. We define the "pre-industrial" baseline as 1850-1865 rather than 1750 (as in the IPCC AR5). We therefore may underestimate the forcing in regions that industrialized relatively early, such as northern Europe. Our focus only on sulfate mass may miss contributions from other aerosol types, like organics. As noted above, the CMIP5 models analyzed in Zelinka et al. (2014) found a substantial contribution to the total ERF ACI from non-sulfate aerosols. However, there is evidence that non-sulfate aerosol levels in the pre-industrial era are underestimated in climate models, and therefore the present-day increase would be overestimated (Gordon et al., 2016).
Additionally, we decided to group the somewhat ambiguous "cumulus rising under stratocumulus" classification with the cumuliform rather than stratiform types, mainly because one of our leading hypotheses for the difference in response in these regions is related to boundary layer decoupling. If this classification were instead treated as stratiform type, our estimates would be larger magnitude overall.
If it turns out that cloud fraction adjustments are negligible in the southeast Atlantic during austral spring due to the persistently overcast baseline but indeed are important in other regions and/or seasons in which the cloud fraction effect would not be saturated , we would expect ERF ACI values to be affected outside the core stratocumulus decks. (If they exist elsewhere, we would expect the cloud adjustments to enhance cloud fraction, although we cannot rule out small reductions in cloud cover.) Finally, our "diurnal" averages using daytime Terra and Aqua data are more heavily skewed toward the morning and early afternoon over the late afternoon and evening. Based on our liquid water path results, we would expect this to lead us to underestimate the countervailing effect of liquid path decreases on Twomey-style brightening. However, averaging together the cloud albedo change estimates from the Terra/Aqua-subset CERES SYN-MHour data actually underestimates the cloud albedo change from CERES EBAF, which should be representative of the daily mean response. It is thus not entirely clear in which direction the bias from not representing a full diurnal cycle should lie.

Conclusions and Future Directions
We find the first observational evidence for substantial cloud microphysical and radiative impacts from shipping on climate-relevant spatial and temporal scales by using universal kriging to estimate cloud properties under the counterfactual situation of no shipping effects in the southeast Atlantic basin. Total albedo changes in the stratocumulus region are dominated by changes in cloud albedo, not in fractional cloudiness. This cloud brightening is driven by an increase in droplet number and decrease in droplet size and is partially offset by a decrease in liquid water path (at least in the afternoon).
Our results join a growing literature of studies using "natural experiments" like volcanic eruptions and the inadvertent modification of clouds from international shipping, which overall have tended to find that climate model estimates of the effective radiative forcing due to aerosol-cloud interactions are likely too negative due to overly strong adjustments that reinforce the Twomey effect (Gryspeerdt, Goren et al., 2019;Malavelle et al., 2017;Toll et al., 2017Toll et al., , 2019. These previous papers, however, have generally taken a "bottom-up" approach of aggregating statistics of individual pollution tracks to reach their conclusions, whereas we take a "top-down" approach to study aerosol-cloud interactions using climatological fields. However, even though these "bottom-up" studies may include hundreds or thousands of individual pollution tracks, they are very likely missing a substantial fraction of pollution-affected clouds due to detectability issues, and it is not obvious that clouds in clearly detectible tracks respond in the same way as clouds in situations in which detection is more difficult (Possner et al., 2018). As automated methods of identifying individual ship tracks improve (Yuan et al., 2019) and other advances are made in tracking air masses affected by individual ships (Gryspeerdt, Smith, et al., 2019), these "bottom-up" and "top-down" methods should be able to converge and provide increased certainty about the net effect of cloud adjustments to aerosol perturbations.
Although the results from the "natural experiments" studies have tended to find small or countervailing cloud adjustments to the Twomey effect, other recent work attempting to control for aerosol-meteorology correlations by stratifying observations by cloud geometrical thickness has come to the opposite conclusion: That cloud adjustments to the Twomey effect are strong and tend to reinforce cooling (mainly due to increased cloud fraction), so much so that compensating warming effects of aerosols on deep convection must be invoked to maintain a realistic global energy balance (Rosenfeld et al., 2019). Based on the relative susceptibilities of cloud fraction and liquid water path (λ) from (the revised version of) Table 1 of Rosenfeld et al. (2019) for Δθ > 18 K, we would have expected a cloud fraction increase of 1.7% (95% confidence interval: 0.9 to 2.6%) and a liquid water path change of 0.2% (−0.7 to 1.0%) for the increase in cloud droplet number concentration from MODIS/Aqua in the subtropical region. (95% confidence intervals are calculated by multiplying the 5,000 Ship-NoShip relative N d change estimates by random draws from a normal distribution with the mean and standard deviation of the susceptibility values.) Both of the estimates based on Rosenfeld et al. (2019) are inconsistent with the 95% confidence intervals of the MODIS/Aqua daytime cloud fraction and liquid water path estimates we found. If the true cloud fraction response in the subtropical southeast Atlantic were as large as suggested by the Rosenfeld et al. (2019) analysis, it seems unlikely that our method would not have been able to detect it. Whether other regions and/or seasons with lower baseline cloud fractions respond in the same manner as those in the southeast Atlantic in austral spring remains an open question, however. Overall, the "natural experiment" literature has thus far provided less insight into changes in fractional cloudiness than in liquid water path-a more dedicated analysis of cloud fraction changes in response to perturbations like shipping and volcanoes could be a promising avenue for future work.
Extrapolating our results globally, we calculate an effective radiative forcing due to aerosol-cloud interactions in low clouds of −1.0 W/m 2 (−1.6 to −0.4 W/m 2 ). Cloud adjustments that decrease liquid water path during the day may reduce the ERF ACI from stratiform-type clouds by a third or more compared to the radiative forcing from the Twomey effect alone.
The uncertainty in our global estimates is dominated by cumuliform clouds even though the central estimates of the forcings are about twice as large for stratiform clouds. Given the plausibility of either strong increases (Albrecht, 1989) or decreases (Seifert et al., 2015) in liquid water path in these clouds, the apparent difficulty of discerning observational signals in these clouds (Grosvenor et al., 2018;Peters et al., 2014), and the general absence of clear pollution tracks within this regime (Toll et al., 2019), better understanding cumuliform cloud responses should be a priority for those interested in better constraining ERF ACI .
The unique setup in the southeast Atlantic provides an opportunity to study climate-scale cloud and radiative responses to a fairly well-known aerosol perturbation in which the causal relationship between aerosol and cloud changes is clear. Cloud microphysical changes in response to shipping emissions and the resulting cloud adjustments in this region could provide valuable benchmarks that can be used to evaluate the representation of the Twomey effect and cloud adjustments in regional and global climate models and in competing observational strategies.