Effect of contrail overlap on radiative impact attributable to aviation contrails

Condensation trails (“contrails”) which form behind aircraft are estimated to cause on the order of 50% of the total climate impact of aviation, matching the total impact of all accumulated aviationattributable CO2. The climate impacts of these contrails are highly uncertain, in part due to the poorly20 understood effect of overlap between contrails and other cloud layers. With the airline industry projected to grow by approximately 4.5% each year over the next 20 years, instances of contrail overlap are expected to increase, including any potential mitigating or amplifying effects on contrail-attributable radiative forcing. However, the impacts of cloud-contrail overlaps are not well understood, and the effect of contrail-contrail overlap has never been quantified. In this study we develop and apply a new model of 25 contrail radiative forcing which explicitly accounts for overlap between cloud layers. Cloud-contrail overlap is found to be responsible for 93% of net radiative forcing attributable to 2015 contrails. We also find significant variation in the sensitivity of contrail radiative forcing to cloud cover with respect to geographic location. Clouds significantly increase warming at high latitudes and over sea, transforming cooling contrails into warming ones in the North-Atlantic corridor. Based on the same data, our results 30 indicate that disregarding overlap between a given pair of contrail layers can result in longwave and shortwave radiative forcing being overestimated by up to 16% and 25% respectively, with the highest bias observed at high optical depths (> 0.4) and high solar zenith angles (> 75°). When applied to https://doi.org/10.5194/acp-2020-181 Preprint. Discussion started: 5 May 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Condensation trails ("contrails") are ice clouds which form in aircraft engine exhaust plumes. Contrails cause both cooling effects through reflecting incoming shortwave solar radiation, and warming effects through absorbing and re-emitting outgoing terrestrial radiation. Previous studies have found the latter 45 effect to be dominant, particularly at night when the cooling effects associated with reductions in incoming shortwave radiation do not exist (Liou, 1986;Meerkötter et al., 1999). The difference between these two effects is the contrail-attributable radiative forcing, which is responsible for long-term climate effects (Penner et al., 1999;IPCC 2013).  Olsen et al., 2013;2: From Ponater et al., 2002, which reports on the same data, 3: Maximum-random overlap is defined in Geleyn and Hollingsworth (1978) by assuming clouds in adjacent layers maximally overlapping and clouds separated by one or more clear layer randomly overlapping)

Source
Target year The net radiative forcing impacts of contrails have been quantified using both global climate models 55 (Chen and Gettelman, 2013;Ponater et al., 2002) and dedicated modeling approaches such as the Contrail Cirrus Prediction Tool (CoCiP) (Schumann, 2012) and the Contrail Evolution and Radiation Model (CERM) (Caiazzo et al., 2017). These approaches have resulted in estimates of total contrail radiative forcing ranging from +15.2 mW/m 2 (Chen and Gettelman, 2013) to +63.0 mW/m 2 (Schumann et al., 2015) for 2006, as shown in Table 1. Normalizing by the total aviation fuel burn in each given year, this and Hollingsworth, 1978). Fromming et al. (2011) also implemented contrails as an increase in cell cloudiness, but assumed simple random overlap with existing clouds. The same assumption is used in 90 Rap et al. (2009), with the HadGEM2 climate model. Burkhardt and Kärcher (2011) assume maximumrandom overlap, as do Burkhardt (2016, 2019). Chen and Gettelman (2013) also implemented contrails in the CAM5 model, representing contrails as an increase in the 3-D cloud fraction. However, they assumed zero overlap between contrails when calculating the impact of linear contrails. These multiple representations are partly responsible for the uncertainty in contrail radiative impact shown in 95 Table 1. However, none of these studies quantified the effect of contrail-contrail overlap.
An analysis of year-2015 global contrail coverage simulated at a resolution of 0.25°×0.3125° using the CERM modeling tool (Caiazzo et al., 2017) provides an estimate of overlap frequency. Assuming maximum overlap by area (i.e. all contrails in a given column overlap to the greatest possible extent, see 100 Section 2.2.2), up to 15% of all contrail area includes overlap with other contrails (Fig. 1, lower plot).
More details on this assumption and the CERM modeling tool are given in Section 2.2. The majority of this overlap occurs for contrails which are no longer line-shaped, and which may appear to be natural cirrus when viewed from the ground. If we exclude contrails which are more than an hour old or which are "sub-visible", having an optical depth below 0.03 (Kärcher, 2002;Kärcher, 2018), this fraction falls 105 to 2.2%. , assuming no contrail-contrail overlap. Lower panel: yearly average coverage (in %), assuming "maximum overlap" such that all contrails in a single column are centered in each 0.25°0.3125° grid cell (in %). Contrail data were generated using the CERM global contrail modeling tool 110 (Caiazzo et al., 2017), which provides contrail quantities and properties discretized to the aforementioned global grid. More information on CERM can be found in Section 2.2.1. Maximum contrail overlap assumes that all contrails in a single vertical grid column overlap to the greatest possible extent by area. This estimate includes contrails which are diffuse and/or "sub-visible" (optical depth < 0.03).
Additionally, there is little agreement on even the sign of how cloud-contrail overlaps change net RF 115 impacts of contrails, due to uncertainty over whether they would more strongly mitigate shortwave (cooling) or longwave (warming) forcing. While Ponater et al. (2002) and Radel and Shine (2008) found evidence for cloud-contrail overlap to reduce the net RF of contrails by 8 to 43%, results by Minnis et al. (1999) and Schumann et al. (2012) implied that clouds could have the opposite effect. This uncertainty is important in light of future changes in cloud cover due to climate change, and projected changes in 120 global patterns of aviation traffic.
Section 2.2.2 includes a broad evaluation of contrail-contrail and cloud-contrail overlap assumptions in literature, as well a detailed description of the assumptions used for this work. https://doi.org/10.5194/acp-2020-181 Preprint. Discussion started: 5 May 2020 c Author(s) 2020. CC BY 4.0 License. In this study, we develop and apply a cloud radiative forcing model designed to address multiple cloud layer overlap in the context of contrails. We first perform a parametric analysis for a single case of overlap between cloud layers, quantifying the potential mitigating or amplifying radiative forcing effects of cloudcontrail and contrail-contrail overlaps. We also estimate the bias in existing estimates which neglect contrail-contrail interactions. This single-column analysis is extended to a global analysis of how contrail 130 radiative forcing is affected by overlap, including sensitivity to season, location, and local conditions (e.g. temperature and solar zenith angle). Finally, we apply this model to estimate the global RF impacts of contrails when considering both cloud-contrail and contrail-contrail overlaps in 2015.

Method
The modeling approach is based on a radiative transfer model previously developed to simulate natural 135 clouds, which we extend to simulate multiple contrail cloud layers. Section 2.1 describes the model and validates the results against existing approaches, and Section 2.2 describes the input data. Using this model, we develop a series of simulations -described in Section 2.3 -which quantify the net radiative forcing impacts of contrail-contrail overlaps and cloud-contrail overlaps under different conditions.

The radiative forcing model 140
The net radiative forcing (RF) attributable to contrails is the sum of two components: longwave (LW) and shortwave (SW). Shortwave radiation is the incoming radiation flux from the sun, which typically undergoes scattering and reflection with minimal atmospheric absorption or re-emission. Longwave ("terrestrial") radiation is the emission of longer-wavelength radiation by the Earth, which undergoes minimal scattering or reflection but is strongly absorbed by clouds before being re-emitted. Contrail cloud 145 layers induce a negative shortwave RF during the day since they reflect incoming solar radiation, slightly increasing the global mean albedo. However, as in the case of natural cirrus clouds, the longwave RF impacts of contrails during both day and night are positive. This is because they absorb terrestrial radiation and re-emit it at the lower temperatures of the upper troposphere (Penner et al., 1999).
In this study we extend and use a cloud radiative transfer model first described by Corti and Peter (2009) which can be applied to both natural or artificial cloud layers (e.g. contrails). This model calculates the cloud-induced change in outgoing longwave and shortwave radiation based on simulated or observed surface conditions (albedo and surface temperature), outgoing longwave flux, meteorological data (ambient temperature), and cloud coverage. The radiative forcing (RF) attributable to a single cloud layer 155 is calculated by performing two simulations: one with the cloud layer present, and one without. The instantaneous RF of a cloud layer is then defined as the difference between the net radiative flux at the top of the atmosphere with and without the layer (IPCC, 2013), so a positive net radiative forcing impact implies an increase in the net energy of the Earth-atmosphere system. 160

Summary of the single cloud layer RF model
The radiative forcing model quantifies the instantaneous RF per unit area of cloud layer. A full description is given in the original model description paper (Corti and Peter, 2009), but we give a brief summary here.
In the original model, the longwave RF is calculated in W/m 2 for a single cloud layer as 165 where L is the outgoing longwave radiation from the surface of the Earth (in W/m 2 ); Lc is the total outgoing longwave radiation from the cloud (in W/m 2 ); Tsrf is the temperature of the Earth's surface (in K); Tc is the cloud temperature (in K); ε is the contrail emissivity; and σ * (the adjusted Stefan-Boltzmann 170 constant, in W/m 2 /K -2.528 ) and k * (= 2.528) are constants (Corti and Peter, 2009). Therefore * * represents the longwave radiation emitted by the cloud (in W/m 2 ) accounting for CO2 and water vapor absorption from the atmosphere (Corti and Peter, 2009). The double-counting of atmospheric absorption is inherent to the original model. More information can be found in Section 3.3.4. Additionally, due to limitations of this approach identified by Lolli et al. (2017), we use an alternative approach for estimating 175 the outgoing longwave flux at the Earth's surface (Section 2.1.2). https://doi.org/10.5194/acp-2020-181 Preprint. Discussion started: 5 May 2020 c Author(s) 2020. CC BY 4.0 License.
The shortwave RF is calculated as where S is the incident solar radiation (in W/m 2 ); is the cloud reflectance for direct radiation; ′ is the cloud reflectance for diffuse radiation; ⍺ is the albedo of the Earth; and t is the atmospheric transmittance above the cloud, assumed constant at a value of 0.73 (Corti and Peter, 2009). The adjusted constants (σ * and k * ) and daily mean atmospheric transmittance (t) are based on clear sky simulations combining results from a high-fidelity radiative transfer model and ECMWF ERA-40 atmospheric profiles (Fu and Liou, 185 1993;Corti and Peter, 2009). Assuming a constant transmittance may result in some bias, as the parameter t would likely vary with location, time and atmospheric composition, including column concentrations of water vapor and aerosols (Schwarz et al. 2020).
While most of the parameters previously mentioned describe the atmospheric conditions, three parameters 190 describe the interaction between cloud and radiation: longwave emissivity (ε) and shortwave reflectances ( and ′ ). All three are dependent on the layer optical depth . Shortwave reflectances, representing cloud interaction with sunlight, are additionally dependent on cloud layer microphysics through the asymmetry parameter g and ′ is additionally dependent on the solar zenith angle. A full description of this derivation is given in Corti and Peter (2009). 195 The optical properties of contrail ice crystals are represented in the model by the asymmetry parameter g of the layer. g measures the degree of anisotropy of scattering and is dependent on the radius and shape of the particle mixture. It ranges from -1 (total backscatter) to +1 (total forward scatter), while equaling 0 for perfect isotropic scattering (Stephens et al., 1990). Ice cloud particles have complex scattering phase 200 functions (Liou et al., 1998;Baran, 2012) but typically fall into the Mie scattering regime with a dominant forward scattering peak, corresponding to an asymmetry parameter between 0.7 and 0.9 (Baran, 2012;Nousiainen and McFarquhar, 2004;Yang et al., 2003). The effect of uncertainty in the asymmetry https://doi.org/10.5194/acp-2020-181 Preprint. Discussion started: 5 May 2020 c Author(s) 2020. CC BY 4.0 License. parameter on contrail RF is investigated in a complementary study (Sanz-Morère et al., 2020). We here assume an average contrail asymmetry parameter, based on in-situ measurements, of 0.77 with a slight 205 increase for the first hour (g = 0.78) (Febvre et al., 2009;Gayet et al., 2012;Schumann et al., 2017). For natural clouds, the asymmetry parameter is calculated as a function of altitude only. We assume that clouds below 8 km have an asymmetry parameter of 0.85 (typical of liquid water clouds); that clouds above 10 km have an asymmetry parameter of 0.7 (typical of long-loved cold cirrus clouds); and that clouds between 8 and 9 km have an asymmetry parameter of 0.8 (Gerber, 2000;210 Kokhanovsky, 2004;Schumann et al., 2017)

Modification, limitations, and validation of the radiative transfer model
We have modified the original approach described by Corti and Peter to account for limitations highlighted by Lolli et al. (2017). Firstly, the original model estimates outgoing longwave flux at the surface by applying a fixed relationship between surface temperature and emitted radiation, based on data 215 from the ECMWF ERA-40 meteorological product. Lolli et al. (2017) found that, below surface temperatures of 288 K, this yielded results that agreed (within 6%) with those from the more complex Fu-Liou-Gu radiative transfer model (Liou, 1986;Fu et al., 1997). However, they also found that for surface temperatures greater than 288 K, this approach is inaccurate and results in radiative forcing errors of approximately 65%. 220 To overcome this issue, we instead use a "top-down" approach in which radiative forcing (longwave) is calculated as the difference between the estimated top-of-atmosphere longwave flux under "clear sky" conditions (without clouds), and the longwave flux perturbed by cloud layer(s). The estimated outgoing longwave radiation in the absence of clouds ( ) is provided in the CERES data product (see 225 Section 2.2.3). This value of outgoing longwave radiation is estimated to have an error of approximately 1.7% (Loeb et al., 2018). However, we do not propagate this error further through our calculations. Hence, we calculate longwave radiative forcing due to contrails as

230
with all other terms as described in Equation (1).
Shortwave radiative forcing is calculated by assuming a constant atmospheric transmittance above the cloud layer, which may result in inaccuracy when considering clouds at different altitudes. This constant value is calculated based on average estimates from a high-fidelity radiative transfer model (Fu and Liou, 235 1993). Lolli et al. (2017) found that the error due to this assumption was negligible, and so we retain it in our model.
We also assume that cloud layers are sufficiently large that 3-D effects (due to horizontal propagation of radiation) are negligible. The effect of this assumption has been investigated in detail previously (Gounou 240 and Hogan, 2007;Davis and Marshak, 2010;Barker et al., 2012;Hogan and Shonk, 2013). Due to the low thickness of contrails, the resulting error in RFSW and RFLW is expected to be on the order of 10% (Gounou and Hogan, 2007).
To ensure that our conclusions are realistic, we also validate the model through comparison to two existing 245 radiative transfer model developed for cirrus clouds: the "Fu-Liou" model (hereafter FL96) (Liou, 1986;Fu and Liou, 1993;Fu, 1996) and CoCiP . We calculate the radiative forcing due to an isolated contrail layer while varying multiple parameters: contrail optical depth, surface albedo, and solar zenith angle (with fixed radiation data). A full description and evaluation is given in Appendix A.
Each of the three models uses a different approach to represent the optical properties of the ice crystals, 250 so initial comparisons are performed by comparing the results when sweeping over a range of input parameters. We find that, for the radiative forcing due to a single contrail, our results match those from CoCiP, with differences of less than 10% for both RFLW and RFSW. Qualitatively, for the same range of particle sizes, FL96 shows similar behavior. However, the magnitude of the calculated radiative forcing differs between our model and FL96, with inconsistencies of up to 40% in RFSW. 255 Due to the strong dependence of RFSW on crystal size and shape, and due to the different treatment of these properties in the three models tested, we conduct a deeper analysis on the resulting difference on RFSW. We choose a specific crystal size in FL96 and compare the simulated RF against results from our model using an "equivalent" asymmetry parameter (more information can be found in Appendix A). For 260 a given surface albedo, we find differences of less than 15% at low solar zenith angles, increasing up to 20% at solar zenith angles greater than 80°. This is consistent with prior evaluations of the two-stream approximation used in our model (Coakley and Chylek, 1975;Corti and Peter, 2009), which has reduced accuracy at high solar zenith angles (see Section 3.3.4). The dependence of RFSW on albedo is also evaluated in each model. Qualitatively the three models show the same behavior with changing albedo, 265 optical depth and solar zenith angle. For albedos below 0.3 the models agree to within 10%, and for albedos below 0.5 the maximum difference is less than 30%. The percentage difference is insensitive to optical depth (see Fig. A2).

Extension to multiple layers
To quantify the effect of cloud-contrail or contrail-contrail overlaps, we extend the model to account for 270 multiple overlapping layers. Computation of longwave RF is accomplished by working outwards from the Earth's surface, as shown in Fig. 2, with each layer absorbing some fraction ε i of the incident longwave radiation while re-emitting a total flux of ε i σ * * . This approach assumes each cloud layer to be at the temperature of the surrounding atmosphere, so that temperature feedbacks can be disregarded and longwave radiation absorption and re-emission is derived from local temperature and surface 275 temperature. Downward fluxes are not shown because the approach neglects temperature feedbacks. As a result, only outgoing radiation is used in our RF calculations.  To calculate the shortwave RF, we start by estimating the shortwave radiation impact of each cloud layer.
Per unit of direct incident shortwave radiation, a fraction of shortwave radiation is reflected and (1 -) is transmitted (absorption of shortwave radiation is assumed to be negligible). The same approach is taken for diffuse shortwave radiation, this time using the parameter ′ . The parameters and ′ are 285 calculated as where is the optical depth of the cloud layer, μ the cosine of the solar zenith angle θ, and γ = 1/(1-g) where g is the layer asymmetry parameter. 290 Due to the high degree of forward scattering of clouds and contrails (Baran, 2012;Nousiainen and McFarquhar, 2004;Yang et al., 2003;Kokhanovsky, 2004), we further assume that (i) shortwave radiation, which has not yet impinged on the Earth's surface, is direct; and (ii) any shortwave radiation reflected from the Earth's surface is diffuse (Corti and Peter, 2009). With these assumptions, the total 295 radiative forcing of two overlapping layers with identical asymmetry parameters is analytically equal to the radiative forcing of a single layer with an optical depth equal to the sum of that from both layers. A full derivation of this result is given in Appendix B (Section B1) for any number of layers.
To model the shortwave radiation impacts of multiple layers, we then collapse the cloud layers into an 300 equivalent single effective layer. To characterize this layer, we derive the effective asymmetry parameter of the overlapping system (Appendix B, Section B2). For N overlapping layers, this is calculated using the optical depth-weighted average value of the gamma function where τi and γi are the optical depth and gamma function (1/(1-gi)) respectively for each individual layer.
Using the effective gamma function, we can then derive and ′ as shown in Eq. (4) and (5) for the full stack of overlapping layers. Substituting (4), (5) and (6) back into Eq.
(2), we obtain the radiative forcing components for N overlapping cloud layers as and then, this can be combined with the previously mentioned procedure for RFLW ( (Caiazzo et al., 2017). CERM follows a bottom-up approach for simulating contrails by 315 combining externally-provided meteorological and atmospheric data with flight track data.
With an hourly time-discretization, and a 0.25°×0.3125°×22 global grid, CERM estimates individual contrail properties (including optical depth and size) for all flights in a year using flight track and atmospheric composition data. CERM models contrails from formation to sublimation based on the 320 physical evolution defined in Schumann (2012). Therefore, it is in theory capable of capturing linear contrails and contrail cirrus. Two contrails allocated in the same grid cell are assumed to be a single contrail while no interaction is assumed between contrails located at different vertical levels. Additionally, physical interactions between simulated contrails and natural clouds are not considered by CERM. This is in part because the contrails may form in the "non-cloudy" parts of grid cells, and in part because of 325 uncertainty over contrail formation when flying through (for example) sub-visible cirrus. We use meteorological reanalysis data from the GEOS forward processing (GEOS-FP) product, supplied by the NASA Global Modeling and Assimilation Office. Flight track data is provided by the Federal Aviation Administration's Aviation Environmental Design Tool.

330
The CERM version used to create the input data for this analysis incorporates new capabilities compared to previous versions (Caiazzo et al., 2017): a higher-resolution vertical grid (22 layers instead of 10 layers); a 4 th order Runge-Kutta advection scheme; and an improved ice crystal coagulation model (Schumann, 2012).

Contrail-contrail and cloud-contrail overlap definition
CERM does not provide the position and orientation of contrails within each grid cell. As such, contrail overlap is computed by assuming the maximum possible overlap, which provides an upper-bound estimate of total overlap. This approach assumes that the smallest contrail (by area) in each vertical column is maximally overlapped with all other contrails in the column, repeating the process for all 340 subsequent contrails in the column. Additionally, CERM aggregates contrails which originate within the same vertical level (at a resolution of ~350 m at cruise altitude) into a single layer. Overlap between contrails which form in such proximity is therefore not explicitly resolved. The same approach is used to model cloud-contrail overlaps, as if clouds were centered in the grid cell.

345
Typically in radiative forcing calculations in literature, clouds and contrails have been assumed to maximally overlap, either by reducing radiation reaching contrail layers or by modeling contrails as an increase in cloud fraction. This is consistent with the fact that cloud coverage is in general larger than contrail coverage (linear in most cases) therefore implying contrail area totally interacting with natural clouds. Contrail-contrail overlaps have been assumed to happen randomly in most of climate models, and 350 maximally in this work and in CoCiP's RF calculations . In this case, additional information exists and can be used to estimate contrail orientation. Due to the existence of flight corridors, overlapping flight paths might be common resulting in potential maximum overlap from contrails for example in the North-Atlantic Corridor. However, this might not happen in denser flight areas like mainland US. Using information on flight paths to include contrail orientation in contrail modeling tools 355 would be useful to more accurately model the impact of contrail-contrail overlaps on contrail-attributable radiative forcing.

Atmospheric radiation data
All atmospheric data required in the radiative forcing model are taken from observations by CERES 360 instruments on three orbiting platforms (NASA Langley Research Center Atmospheric Science Data Center 2015). CERES data are provided at three-hour intervals, so hourly average values are used.
The terrestrial, longwave radiation flux is simulated using the estimated "clear-sky" outgoing longwave flux provided by CERES. The "clear sky" flux is the estimated flux in the absence of clouds. This removes 365 the need to estimate outgoing fluxes based on indirectly-observed surface temperatures. Longwave emission from cloud layers is calculated as described in Corti and Peter (2009). The total incident shortwave radiation S is computed using the solar zenith angle calculated based on time and geographic location (Kalogirou, 2014) as 370 where S0 is the solar constant (1366.1 W/m 2 ), µ is the cosine of the solar zenith angle θ, and J is the Julian Day.

Natural cloud data
CERES instruments also provide data on natural cloud coverage. These are divided into four vertical levels defined by pressure, and include cloud properties such as optical depth and temperature. We use 375 this data to estimate natural cloud cover when calculating the impacts of contrails in 2015.
Since CERES instruments observe both contrails and natural cirrus clouds, we may be double counting the influence of contrails. Four levels of clouds are given in CERES data, defined by their pressure level and corresponding to the following altitudes: from 0 to 10,000 ft, from 10,000 ft to 16,500 ft, from 16,500 380 ft to 30,000 ft, and above 30,000 ft. Accordingly, most contrails would appear in the 4 th level detection.
There is a high-level cloud in the same location as a CERES detectable contrail (optical depth higher than 0.02) in 58% of contrail cases, whereas only 6% of contrails occupy the same level as mid-level clouds (3 rd CERES vertical level). There is in theory the possibility that ~60% of all contrails are already 385 accounted for in the CERES data. However, the detection limit of the CERES instruments is approximately τ = 0.02 (Dessler and Yang, 2003). Considering that the average optical depth from CERM for 2015 global contrails is 0.065, a significant fraction of CERM contrails are not detectable by CERES instruments.

Experimental design 390
We analyze the radiative forcing impacts of cloud-contrail and contrail-contrail overlaps using a threestep approach.
In the first step, through a parameterized analysis, we quantify the effect of a two-layer overlap on total radiative forcing when compared to a case where the layers are assumed to be independent, calculating 395 how the effect of overlap varies as a function of the layer properties and the local conditions. This analysis shows the conditions under which the RF of two overlapping contrails is significantly different to the total RF of two independent contrails. This nonlinearity is representative of the existing bias in estimated contrail RF values found in the literature.

400
In the second step, we evaluate the global sensitivity of contrail RF to cloud-contrail and contrail-contrail overlaps using 2015 atmospheric data (meteorology and natural clouds). We calculate the RF associated with one or two contrail layers at each global location for one day from each month of the year in order to capture seasonal variation. To demonstrate this, we simulate a case used previously in estimates of contrail-attributable radiative forcing (Myhre et al., 2009;Schumann et al., 2012). The RF attributable to 405 a hypothetical contrail is calculated for each location globally assuming typical optical properties (g = 0.77), optical depth (0.3), and altitude (around 10.5 km). In order to quantify the effect of cloud overlap we evaluate radiative forcing with and without natural cloud cover ("all sky" vs. "clear sky"). By subtracting the RF obtained in the "clear sky" scenario from the RF obtained in the "all sky" scenario, we obtain the difference in contrail RF attributable to the presence of clouds. The results can then be linked 410 to different cloudiness conditions to systematically analyze the impact of cloudiness on contrail RF. In order to quantify the global sensitivity to contrail-contrail overlaps we simulate a superposition of two contrail layers at each location, separated by a vertical distance of approximately 0.5 km.
Finally, we quantify the effect of cloud-contrail and contrail-contrail overlap on the global contrail-415 attributable RF in 2015. We use year-2015 contrail coverage data obtained from CERM (Caiazzo et al., 2017) and analyze the associated radiative forcing impacts for the four scenarios shown in Table 2.

Overlapping (O)
OC OA Global evaluations are performed using detailed contrail coverage estimates and meteorological data 420 described in Section 2.2.

The effect of overlap on contrail-attributable radiative forcing in a single column
In this section we evaluate the general effect of overlap on contrail RF through a parameterized analysis.
We simulate two overlapping layers with different optical depths () and temperatures (T) (either natural 425 cloud or contrail). By varying the layer properties, we are able to simulate both cloud-contrail and contrail-contrail overlaps. We also evaluate the effect of solar zenith angle (θ), estimated outgoing longwave radiation without clouds (OLRclear), and Earth surface albedo (α).
The contrail modeling and observation literature suggests that contrails are usually optically thin, with 430 typical optical depths in the range 0.05 to 0.35 (see Table 1). They also form almost exclusively at cruise altitude. Natural clouds are located within a greater range of altitudes and can achieve greater optical depths. We simulate contrail layers over a range of depths (0 <  < 0.5), based on typical values, at low temperatures/high altitudes (210 -230 K), and with an asymmetry parameter of 0.77, representative of mature contrails (Schumann et al., 2017). Cloud layers are simulated as being thicker (0 <  < 4), at higher 435 temperatures/lower altitudes (215 -280 K), and with an asymmetry parameter of 0.85, corresponding to low level clouds. When not otherwise specified, we assume each contrail layer to have an optical depth  of 0.3 and temperatures of 215 K (upper) and 220 K (lower). This optical depth is at the upper bound of literature estimates of typical values for contrails (Voigt et al., 2011). For this analysis natural cloud layers are assumed to have an optical depth  of 3 and a temperature of 260 K. The prescribed outgoing longwave 440 radiation in this single-column analysis is 265 W/m 2 (consistent with a ~288 K surface temperature), with an albedo α = 0.3 and solar zenith angle θ = 45º.
The total forcing for the combined, overlapping layers is calculated as shown in Section 2.1. We calculate the "independent" forcing as the RF that would have been calculated by adding together the RF from each 445 layer independently, without accounting for any overlap. We evaluate the effect that overlap has on the contrail-attributable net radiative forcing in both systems (cloud-contrail and contrail-contrail) as a function of each parameter. We then estimate the bias in estimated RF that results if overlap is ignored, as is typically the case in existing contrail modeling. We also evaluate contrail RF when surrounded by cirrus clouds and finally, we validate our overlap model comparing it with the FL96 model previously 450 described (Liou, 1986;Fu and Liou, 1993;Fu, 1996).

Parametric analysis of cloud-contrail and contrail-contrail overlap effects on contrailattributable net RF
The effect of overlap on contrail-attributable RF depends both on cloud layers' properties and on local conditions. We first evaluate how the effect of overlap varies with cloud layer properties, including 455 thickness of the two layers. We then quantify the effect of local conditions: solar zenith angle (θ), estimated outgoing longwave radiation in clear sky conditions (OLRclear), and Earth surface albedo (α).

465
We evaluate the effect of overlap on contrail-attributable net RF for both cloud-contrail (with the contrail at 215 K) and contrail-contrail (at 215 K and 220 K) systems. The variation in contrail-attributable net RF with optical depth of either layer is shown in Figure 3. A decomposition of the results in terms of longwave and shortwave components can be found in the SI (figures S1 and S2). The panels on the left 470 show the effects of cloud-contrail overlap, while those on the right show the effects of contrail-contrail overlap. The upper row shows the net RF when the layers are considered to be independent, while the bottom row shows the RF when accounting for interactions between the two (i.e. overlap). Each panel shows the net, contrail-attributable RF of the system (i.e. subtracting only any RF which is calculated when no contrails are simulated). 475 The RF attributable to a single contrail (no overlap) as a function of its optical depth is shown in the upper left panel (Fig. 3a). This is because, when overlap is ignored, the contrail-attributable RF of a cloudcontrail system is equal to the RF of the contrail alone. The RF increases from zero to a maximum of ~1.2 W/m 2 as the optical depth increases to ~0.2, after which increasing depth instead results in reduced RF. 480 This is due to the compensation of the increase in absorption by the increase in reflectance with increasing optical depth. The lower left panel (Fig. 3c) then shows how the presence of a cloud layer affects contrailattributable RF as a function of the optical depth of each layer. The presence of a (lower) natural cloud layer can either increase or decrease the contrail-attributable RF depending on the optical depth of the cloud layer. Thin clouds can transform a warming contrail into a cooling one by absorbing part of the 485 longwave radiation that previously reached the contrail. Thick clouds can transform a cooling contrail into a warming one (from a net RF of -0.54 W/m 2 to +4.1 W/m 2 at a contrail optical depth of 0.5) by mitigating the shortwave benefit of the contrail. These results explain the existing uncertainty related to the effect of natural clouds on contrails' radiative impact. If overlap between the layers is ignored (Fig.   3a), these features are not captured. 490 Figure 3d shows the effect of contrail-contrail overlaps on contrail-attributable RF. The effect of each contrail individually can be seen based on the values along the left and lower edges. The lower contrail, due to its higher temperature (less LW absorption), becomes cooling at a lower optical depth of ~0.22 (compared to ~0.45 for the upper contrail). The effects of overlap are similar to the effects obtained when 495 a thin cloud (τ ~ 0.1) is overlapping with a contrail: the net effect of increasing the optical depth of the contrail is to make the system more cooling (Fig. 3d). However, since both layers are artificial (contrails), increasing the optical depth of either layer yields a more negative RF, unlike the case of a thick natural cloud with a thin contrail. This is because the shortwave benefit attributable to contrails increases regardless of which layer is providing the shortwave benefit. This results in a monotonic decrease in 500 warming (increase in cooling) attributable to the net contrail-attributable RF, from +1.2 W/m 2 for a single contrail of optical depth 0.25, to -10 W/m 2 for two contrails both of optical depth 0.5. For comparison, Figure 3b (upper left panel) shows the result when RF is calculated based on the independent combination of each contrail's RF. Independent calculation gives the wrong response by neglecting the screening effect on longwave radiation by the lower contrail. This error is small for low contrail thicknesses, with a 505 maximum difference of -1.0 W/m 2 for a total contrail-contrail system thickness below approximately 0.15. However, for thicker contrail layers, both the sign and magnitude of the net effect can be incorrectly predicted when overlap is neglected. This analysis also confirms the findings of Kärcher and Burkhardt (2013) with regards to the overestimation of contrail RF by prescribing a mean optical depth. As an example, two simulated overlapping contrails of optical depths 0.1 and 0.2 result in ~0.8 W/m 2 of radiative 510 forcing, but two overlapping contrails of optical depth 0.15 result in a forcing of 1.1 W/m 2 ).
The altitude (temperature) of each layer also affects the effect that overlap has on the net contrail-attributable RF. Net attributable RF of a contrail-contrail system decreases as contrail altitude decreases (increasing temperature), due to the increase in the temperature of re-emission. For a cloud-515 contrail system, the contrail-attributable RF is most sensitive to the altitude (temperature) of the natural cloud. The absolute difference varies from +6.1 W/m 2 , for warmer (lower-altitude) clouds, to -12 W/m 2 , for cooler (higher) clouds, assuming an optical depth of 3 for the natural cloud layer (see Fig. S3 in the SI).

520
The radiative forcing attributable to contrails (as well as the effect of overlap) also varies as a function of local conditions, such as the outgoing longwave radiation (related to surface temperature), surface albedo, and solar zenith angle. The greatest contrail-attributable warming occurs for high values of outgoing (terrestrial) longwave radiation, and high surface albedos. This is due to the combination of increased longwave radiative forcing and the reduced shortwave benefit from the contrail. We also find that the net 525 RF of the contrail-contrail system is reduced as the solar zenith angle increases. As θ increases from 0º to 75º, the maximum net RF (at maximum OLRclear and α) decreases from 27 W/m 2 to 8.0 W/m 2 . This effect, driven by changes in the shortwave cooling, is explored in more detail in Appendix B (Section 3), with additional figures in the SI ( Figure S4). The relative effect of overlap on both the warming and cooling components of contrail-attributable RF is, in relative terms, insensitive to outgoing longwave radiation 530 and albedo. Due to the low absolute values of |RFSW| at maximum α and high values of |RFLW| at maximum OLRclear, maximum absolute net RF decrease happens in those areas.

Bias in existing estimates of contrail RF due to layer overlap
We now evaluate the error in both RFSW and RFLW which results from ignoring the effect of contrailcontrail overlap. We use RFO to denote the RF when overlap is treated explicitly and RFI to denote when 535 overlap is ignored ("independent"), in which case the total RF is the sum of the RF from each cloud layer.
The relative change in the estimated RF impact of the system is then where a positive value of D indicates that the assumption of independence results in an overestimate of 540 warming effects (LW) or an underestimate of cooling effects (SW). Equivalently, a positive value means that accounting for overlap results in a decrease in the RF of the system relative to the independent calculation. Figure 4 shows the percentage bias resulting from ignoring overlap when quantifying the RF of a contrail-545 contrail system. This is quantified as a function of each contrail's optical depth and of the local solar zenith angle (θ). In each case, the upper and lower contrail have identical physical properties, as described in Section 3.1.1. We find that accounting for overlap consistently results in a reduced longwave RF for two overlapping contrail layers. This means that, if overlapping contrails are considered as independent, their longwave RF is overestimated by up to 16% (for contrails with optical depth of 0.5). This effect is 550 independent of the solar zenith angle.  For shortwave RF, the error resulting from independent calculation is sensitive to the solar zenith angle. 560 In most cases, the total shortwave ("cooling") RF is smaller in magnitude when correctly accounting for overlap, relative to the independent calculation. This corresponds to an overestimate of the total reflectance if contrails are treated as independent. The magnitude of this error generally increases with contrail optical depth. Near sunrise or sunset (θ ≈ 75°), accounting for overlap reduces the calculated cooling effect by 25% for τ = 0.5. However, we observe a change in the sign of the error at zenith angles 565 below ~25°. At noon (θ = 0°), assuming independent effects results in a slight underestimate of the cooling effect for any optical depth between 0 and 0.5, up to a value of 3.2%. The cause for the change in sign at very low solar zenith angles is investigated in detail in Appendix C. The effect on total net RF depends on the tradeoff between the effects on both RFLW and RFSW. At low 570 solar zenith angles, neglecting contrail-contrail overlaps results in an overestimation of net RF. Due to the changes in sign of the error for shortwave RF, and the fact that the magnitude of each of the two components varies based on different factors, the effect on net RF at high solar zenith angles will depend on factors such as the location, time, and properties of each contrail.

575
In summary, we find that the net radiative forcing due to contrails may include a significant non-linear term due to overlap which is not captured in existing models. For contrails with optical depths of up to 0.5, we find that failing to account for this non-linearity could result in an overestimate of both the longwave warming (up to 16%) and the shortwave cooling (up to 25%). The sign and magnitude of the effect on the system net RF is highly dependent on layers' properties, local conditions, and the solar zenith 580 angle. The total effect of overlapping on a single contrail is therefore dependent on the solar zenith angle (time), temperature (altitude), and geographic location in which the contrail is formed.

Parametric analysis of radiative impact from a contrail located in-between cirrus clouds
We also model the case of a single contrail located between two natural cirrus cloud layers. We simulate 585 a single contrail with the same properties as were used in the previous section (temperature of 215 K, optical depth of 0.3, and asymmetry parameter of 0.77). This is bracketed by two cirrus clouds, 500 m above and below the contrail, with optical depths of up to 1.5 and an asymmetry parameter of 0.75 (Kokhanovsky, 2004). 590 Figure 5 shows how the single contrail RF varies as a function of the optical depth of both natural cirrus clouds and as a function of solar zenith angle. For reference, the estimated RF for the contrail at a solar zenith angle of 45° in the absence of clouds is +27.9 W/m 2 (longwave) and -26.9 W/m 2 (shortwave), resulting in a net forcing of 1.0 W/m 2 . https://doi.org/10.5194/acp-2020-181 Preprint. Discussion started: 5 May 2020 c Author(s) 2020. CC BY 4.0 License. The presence of either cloud layer alone decreases both the longwave and shortwave RF attributable to 600 the contrail, as previously discussed. Except at high solar zenith angles, increasing the optical depth of either cloud layer reduces the net RF of the contrail layer. This is because the contrail's longwave RF falls rapidly, while the shortwave RF is less affected. The contrail's longwave radiative forcing decreases by up to a factor of seven when the surrounding clouds are sufficiently thick (τ = 1.5), while the shortwave radiative forcing is only reduced by a factor of three. However, at high solar zenith angles, this situation 605 is reversed (see Figure B2 in Annex B). This means that the contrail RF instead initially increases with increasing cloud thickness.

Validation of the overlap model
In addition to validating the model for the purposes of simulating a single contrail (Section 2.1.2), we also 610 compare the model's estimates of the effect of two-layer overlap to estimates from an existing radiative https://doi.org/10.5194/acp-2020-181 Preprint. Discussion started: 5 May 2020 c Author(s) 2020. CC BY 4.0 License. transfer model -the previously-described Fu-Liou radiative transfer model (FL96). FL96 uses solid hexagonal columns to represent ice clouds, which have previously been found to be best represented in the Corti and Peter model by assuming an asymmetry parameter g = 0.87 (Corti and Peter, 2009). Figure   6 shows the error resulting from considering overlapping contrails as if they were independent, for both 615 longwave and shortwave components, in both models. Qualitatively, the behavior is consistent between the two models. Both models estimate that the 625 discrepancy in simulated longwave and shortwave RF (comparing the "overlap" to "independent" cases) increases with the increasing optical depth of each cloud layer. We also observe the same reversal of sign in the shortwave error at very low solar zenith angles. FL96 finds that both errors increase more quickly with optical depth than is estimated by our model, finding a maximum error in longwave RF of 25% (17% in our model) and in shortwave RF of 24% (18% in our model). This indicates that our model correctly 630 represents overlapping behavior but might underestimate the effect on both terms. The net RF difference is always lower than 30% and varies with solar zenith angle. At low solar zenith angles we underestimate net RF (both for two independent and overlapping contrails). At θ = 45º we obtain the best agreement, with differences lower than 10% and at θ = 75º we overestimate net RF by up to 30% at an optical depth, for both contrails, of 0.5 (at the upper end of current contrail optical depth estimates). 635

Global sensitivity of cloud-contrail and contrail-contrail overlap to location and season
We next quantify the variation of contrail radiative forcing as a function of geographic location and time of year. This captures the primary drivers in variations regarding the effects of overlap, as identified previously.

640
To obtain these sensitivities, we run a global simulation using 2015 atmospheric data (including radiation and natural cloudiness data as described in Section 2.2) in which we simulate the presence of a contrail layer in each location across the globe. We here assume that, in each grid cell, 1% of the total area is covered by contrail, reproducing an analysis performed by Schumann et al. (2012). We evaluate the effect of both cloud-contrail and contrail-contrail overlaps on contrail-attributable RF. We also calculate the 645 error which would be incurred by treating two overlapping contrails as independent. Figure 7 shows the radiative forcing per unit of additional contrail optical depth at each location, under both "clear sky" and "all sky" conditions (without and with natural clouds respectively, for year-2015 natural cloud cover). The RF varies as a function of latitude, consistent with prior studies (Schumann, 650 2012). The longwave warming (RFLW) is maximized in regions with higher surface temperatures such as the equator. Cooling (negative RFSW) is instead sensitive to surface albedo, being maximized over oceans and minimized over snow-covered or desert regions.

660
By comparing the "all-sky" and "clear-sky" simulation results, we find that the absolute value of both components of radiative forcing is reduced by the presence of clouds. The global mean reduction in shortwave forcing (~83%) exceeds the reduction in longwave forcing (~42%), meaning that cloud overlap causes a more than three times increase in the global, area-weighted average, contrail-attributable net RF, 665 from +27.8 mW/m 2 to +107.1 mW/m 2 per unit of contrail optical depth. These values are consistent with prior studies (e.g. Schumann et al. 2012).
Our assumed asymmetry parameter for each contrail layer (g = 0.77) corresponds to a greater backscatter than is the case in previous studies (Fu and Liou, 1996;Myrhe et al., 2001;Schumann et al., 2012). This 670 explains the low global sensitivity obtained in clear-sky conditions. For comparison, using an asymmetry parameter of g = 0.9 (typical of regular, spherical particles) results in a global mean, clear-sky sensitivity of +144.3 mW/m 2 , reducing cloud-contrail global impact. A deeper analysis of uncertainty related to microphysics and resulting global sensitivity to contrail is the subject of a complementary work (Sanz-Morère et al., 2020). 675 At night this effect reverses, as the reduction in reflected shortwave radiation is lost while the reduction in absorbed longwave radiation remains. The global, area-weighted average nighttime contrailattributable RF is therefore reduced by 42% when accounting for the presence of clouds. However, these effects vary significantly with geographic location. 680 The depth, frequency, and altitude of natural cloud cover all vary as a function of location, resulting in a geographical dependence of the sensitivity of contrail RF with respect to clouds. Thick, low altitude clouds are more common at midlatitudes, while higher, thinner cirrus clouds are more common in the tropics (Warren et al., 1988;Sassen et al., 2008;Marchand et al., 2010). The effect of these clouds on 685 contrail RF is shown in Figure 8. In the tropics (TROP, 30°S -30°N), contrail RF is 1.5 times higher in the presence of clouds. However, in the midlatitudes (MLAT, 30°N -60°N), the thicker, warmer clouds have a greater effect. Overlap with midlatitude clouds increases the net RF attributable to a contrail by more than a factor of six, from 8.7 mW/m 2 to 66 mW/m 2 . This result is consistent with the analysis given in Section 3.1.1, and is due to the high reflectivity of the thick, low-altitude clouds. 690 We also quantify the sensitivity of contrail RF to overlap in four different geographical subregions: area 1, representing the North Atlantic corridor; area 2, which includes parts of Asia; area 3, approximately representing the continental United States; and area 4, approximately representing Europe (see Figure 8).
These areas include ~53% of all passenger traffic in 2017 (Boeing, 2018) and differences in sensitivity 700 for each region provide insights into the effects of future growth.
In all four regions, clouds have a greater relative and absolute effect on shortwave RF than on longwave RF (Figure 8). In area 3, clouds reduce the longwave RF per unit contrail optical depth by 46%, while reducing the shortwave RF by 83%. This results in a 2.3 times increase in the net RF relative to the clear-705 sky case. By contrast, in the North Atlantic corridor (area 1), clouds reduce the longwave RF by 44%, but the shortwave RF is reduced by 99%. This changes a cooling effect of 70 mW/m 2 into a warming of 690 mW/m 2 . The effects of cloud overlap in areas 2 and 4 lie in between these two extremes.
These variations are driven by differences in natural cloud coverage (primarily due to latitude) and surface 710 albedo (e.g. land vs. sea). In the case of area 1, contrails are mostly forming over water, which has a very low albedo. As a result, there is a larger shortwave benefit, and therefore a greater increase in the net RF when this benefit is removed by overlap with clouds. By contrast, over area 3 there is a greater land fraction and the clouds are thinner, resulting in a smaller overlap effect. These results suggest that avoiding overlap of contrails with clouds will yield the greatest benefit on midlatitude, oceanic routes, 715 whereas the benefits of doing so over land and/or at lower latitudes will be smaller.
Contrail RF, and its sensitivity to clouds, also varies by season. Under all-sky conditions, in the Northern Hemisphere, the net contrail sensitivity is globally 15% lower in local winter than in local summer. This 720 is because the longwave benefit due to cooler surface temperatures exceeds the shortwave disbenefit from shorter days (less insolation). However, this varies significantly by latitude because of the effect of changes in day length.
Climate change is likely to affect these results due to its effects on global cloud cover (Norris et al., 2016). 725 Current satellite data show that cloud top heights are gradually increasing, which will likely decrease contrail net RF due to the resulting decrease in cloud top temperature. It is also anticipated that the tropics will expand (Kim et al., 2017). This will mean that more contrails are overlapping with high-altitude clouds, resulting in a reduced sensitivity to cloud overlap as discussed earlier. 730 We also evaluate how the effect of contrail-contrail overlap on contrail-attributable RF varies by location. This is quantified by simulating two contrail layers at each location, first treating them as independent and then calculating the total RF when accounting for overlap. The layers are simulated as being separated by 500 m. We find that correctly accounting for overlap results in a decrease in both the cooling and warming effects, relative to the "independent" calculation. The percentage decrease in each component 735 is approximately uniform across all locations (consistent with Section 3.1.1). Since the components are of opposite sign, this results in a non-uniform effect on total net RF. Contrail overlap has the greatest effect on the net RF when contrails are located in hot, equatorial areas (increased longwave disbenefit) with high albedo (reduced shortwave benefit), as is the case in low-latitude desert areas such as the Sahara.
This results in a maximum contrail attributable net RF reduction by contrail-contrail overlapping in the 740 tropics (TROP), where we find a reduction from an average sensitivity of 1.6 W/m 2 (per unit of optical depth) for two "independent layers" to an average sensitivity of 0.6 W/m 2 for two "overlapping layers".
Global sensitivity maps to contrail-contrail overlap are shown in Figure S5 of the SI.

Effect of cloud-contrail and contrail-contrail overlaps on net 2015 global radiative forcing attributable to contrails 745
Finally, we quantify the net effect of cloud-contrail and contrail-contrail overlap for existing aircraft traffic patterns. We use year-2015 contrail coverage data as estimated using CERM (see Section 2.2.1).
The RF impacts of contrails are presented in Table 3, under all-sky and clear-sky conditions, and with and without explicit treatment of contrail-contrail overlap.

Cloud-contrail overlaps 750
For 2015, we find that approximately 75% (by area) of contrails overlap with mid-level clouds. We compare results calculated under all-sky and clear-sky conditions (scenarios OA and OC) to quantify the effect of cloud-contrail overlap on contrail-attributable RF. https://doi.org/10.5194/acp-2020-181 Preprint. Discussion started: 5 May 2020 c Author(s) 2020. CC BY 4.0 License. Figure 9 shows the effect of cloud-contrail interactions on the shortwave and longwave radiative forcing 755 due to contrails. We find a 66% decrease in net global cooling attributable to contrails as a result of cloud cover, accompanied by a 37% decrease in warming. Accounting for cloud interactions therefore results in more than ten times greater net contrail-attributable warming. As a consequence, the annual average, global net RF changes from +0.7 mW/m 2 under clear sky conditions to +9.7 mW/m 2 when including clouds ("all-sky"). Clouds are then responsible of 93% of 2015 contrail climate impact. At night, contrails 760 over natural clouds have a lower net RF due to the lack of any shortwave effect. As a result, the presence of natural clouds during nighttime reduces the net RF of contrails by 37% as the only effect that clouds can have at this time is to mitigate the contrail-attributable longwave RF.

Contrail-contrail overlaps
Under an upper-bound assumption for the total area of contrail overlaps, we find that 15% of all modeled contrail area overlaps with other contrails at different altitudes. If the effect of cloud-contrail overlap is 770 ignored, the maximum contrail-contrail overlap results in a more than three times increase in the contrailattributable net radiative forcing. This is made up of a 21% reduction in longwave warming but a 38% decrease in shortwave cooling. However, if cloud-contrail overlap is accounted for, the net impact of contrail-contrail overlap is instead a 3.0% reduction in net contrail-attributable RF. The reduction in longwave warming is 2.0%, exceeding the 1.8% reduction in shortwave forcing. This difference is due to 775 the strong mitigation of shortwave forcing (approximately 1/3 of that under clear sky conditions) by existing clouds, and is consistent with the global sensitivity to contrail-contrail overlaps demonstrated in Section 3.2. The majority of contrail-contrail overlap occurs in low-albedo areas such as the North-Atlantic corridor (area 1) or at high latitudes (areas 3 and 4), resulting in a small absolute effect on net RF (-0.3 mW/m 2 ). 780 These results are sensitive to the assumptions regarding the degree of overlap in each model column. We assume that all contrails in a given model column overlap to the maximum extent, providing an upper bound for the total effect of contrail overlap. If we instead assume minimum overlapwhere each contrail in the column "avoids" overlap until there is no remaining uncovered areathen contrail-contrail overlap 785 only occurs for 2% of all modeled contrail area. This limitation is explored further in Section 3.3.4.  Table 3 shows the total contrail-attributable RF with and without clouds, and either accounting for or 790 neglecting the effects of contrail-contrail overlap. We find that contrails induce a net RF of 9.7 mW/m 2 for 2015. This result includes a 3% reduction in overall RF from contrail-contrail overlap, but most of it (93%) is due to overlap with clouds.

Overall impact of cloud-contrail and contrail-contrail overlap on global RF
These results suggest that the impacts of cloud-contrail overlap are significant, but that contrail-contrail 795 overlap can likely be neglected in modeling studies under current conditions. However, our result of +9.7 mW/m 2 for the net impact of contrails is at the low end of existing literature estimates (see Table 1). This is due to uncertainties in contrail coverage, contrail optical depth, and contrail optical properties. The global CERM simulation output has an average optical depth per contrail of 0.065 and a global coverage of 0.39% by area, both of which are at the lower end of literature estimates (see Table 1). As a sensitivity 800 test, if we increase the optical depth of all contrails from the CERM output data by a factor of four to give the same average per-contrail optical depth as Schumann et al. (2013), which found a net RF of 49.2 mW/m 2 , we find a global net contrail RF of 32.6 mW/m 2 . Under these conditions, we find that contrailcontrail overlaps decrease the simulated global RF by 8%.

Radiative transfer model
Our radiative forcing model is based on an existing, single layer cloud model (Corti and Peter, 2009). This model has some limitations. 810 When calculating the total outgoing longwave radiation for each layer, the model includes an estimate of absorption by atmospheric CO2 and water vapor. Estimates for multiple overlapping layers may therefore double-count this contribution. Additionally, cloud emissivity is estimated as only a function of the cloud optical depth. These limitations may partially explain some of the differences in the calculated outgoing 815 longwave radiative forcing between this model and the Fu-Liou radiative transfer model, as discussed in Section 3.1.4.
The longwave radiative forcing model also assumes all layers to be in equilibrium, and does not account for local temperature feedbacks due to the presence of artificial cloud layers. Finally, we do not account 820 for 3-D effects. Cloud layers are assumed to be vertically homogeneous and edge effects are ignored, as in the reference model. A previous investigation of contrail radiative forcing found that 3-D effects could change simulated radiative forcing by ~10% (Gounou and Hogan, 2007).
Regarding shortwave radiative forcing, we do not account for inhomogeneity in the above-cloud 825 atmospheric transmittance of shortwave radiation, instead considering it to be constant at 73%. Shortwave radiative interactions between contrails and other constituents (such as tropospheric aerosols and water vapor) are also not explicitly accounted for. Secondly, the adapted model uses a two-stream approximation of radiative transfer (Coakley and Chylek, 1975). This has been shown to give accurate results at optical depths below ~1 and solar zenith angles below 75°. However, this results in inaccuracy 830 outside of this range, as shown by comparison to other models (Appendix A). Annually, the solar zenith angle is between 75 and 90° for 16% of the time globally, and 14.5% of the time at latitudes covering the majority of current commercial flights (30 o N -60 o N).
The two-stream approximation used in this model is most accurate for low optical depths. This is 835 appropriate for contrails and thin natural cirrus, but lower-altitude natural clouds can be much thicker.
For this reason, we use an asymmetry parameter for high altitude clouds and contrails based on direct observations (Sanz-Morère et al., 2020), while using an asymmetry parameter similar to that suggested by Corti and Peter (2019) for low altitude clouds.

840
In order to improve the presented radiative forcing model, future research with this approach may wish to prioritize a more accurate model of natural clouds and atmospheric interactions, such as by including more direct satellite observations. 845

Input data
Due to the lack of additional input information, and to provide a conservative estimate, we assume that all contrails overlap maximally within a column. This assumption would not be necessary if additional information was supplied by the base contrail model. Additionally, contrail coverage could be assessed, or confirmed, by satellite measurements. Some studies (Kärcher 2009;Iwabuchi et al., 2012) have 850 combined satellite imagery (e.g. from MODIS) with observed cloud coverage data in order to provide an improved estimate of contrail coverage. The combination of these data with single-contrail modeling tools (such as CERM) may help to improve the accuracy of estimated contrail coverage. However, there remain significant uncertainties due to the non-detection of very thin contrails (Kärcher et al., 2002), as well as the difficulty of distinguishing between between long-lived contrails and natural cirrus clouds in 855 observational data. Finally, the natural cloud data provided by CERES is coarsely resolved into only four layers in the vertical dimension, and lacks some additional useful information. Alternatives to CERES like CALIPSO or CloudSat (Iwabuchi et al., 2012;Tesche et al., 2016) may provide a useful alternative, as they include both more precise estimates of cloud altitude and additional optical properties of the cloud layers. 860 These results are also sensitive to the optical depth of the simulated layers. Contrails simulated by CERM have a mean contrail optical depth of 0.065, at the lower end of a significant uncertainty range based on existing literature (see Table 1). Since the effects of overlap increase non-linearly with optical depth, estimates based on models which predict thicker contrails may find a significantly greater impact of 865 overlap. Finally, there remain significant uncertainties in contrail coverage. Improved estimates of contrail lifetime and formation frequency could significantly affect the frequency, and therefore total impact on contrail-related RF, of cloud-contrail and contrail-contrail overlap.

Conclusions
We develop and apply a radiative transfer model to estimate the effect of cloud-contrail and contrail-870 contrail overlap on the net radiative forcing from contrails. The results will not only improve our understanding of the nonlinearities in global RF impacts and help quantify current and future RF impacts https://doi.org/10.5194/acp-2020-181 Preprint. Discussion started: 5 May 2020 c Author(s) 2020. CC BY 4.0 License. more accurately, but will also help to inform policymakers and researchers to identify technical, operational, and regulatory means to reduce these impacts.

875
We find that overlap between contrails and natural cloud layers can cause a non-linearity in the net radiative forcing. In most cases, overlap between a contrail and a second cloud layer reduces both the cooling and warming effects of the contrail. This effect is sensitive to the optical depth of each cloud layer. We find a net increase in radiative forcing when contrails overlap with thick clouds (τ > 0.5), but a net decrease when contrails overlap with thinner clouds. However, overlap between two contrails is in 880 general beneficial for climate, decreasing the total contrail-attributable RF. The magnitude of this effect is sensitive to local conditions, including surface albedo, solar zenith angle, and surface temperature.
Under night-time conditions, overlapping between contrails and any other cloud layer consistently reduces the net contrail RF due to the lack of competing shortwave effects.

885
The radiative forcing attributable to a contrail layer increases by a factor of three due to the presence of natural clouds on a global mean basis, but this varies by region. Clouds have a greater effect on midlatitude contrail radiative effects than in the tropics due to their greater thickness and lower altitude. They also have greater effects over oceanic routes. We find that contrails over the North Atlantic corridor have, on average, a small cooling effect under clear-sky conditions (-0.07 W/m 2 per unit of optical depth) but cause 890 warming (+0.69 W/m 2 per unit of optical depth) in cloudy conditions, suggesting that avoiding cloudcontrail overlaps in this region could yield climate benefits. This sensitivity also varies by season, with a 15% decrease in RF per unit of optical depth in the Northern Hemisphere from summer to winter.
For year-2015 atmospheric data and flight activity, we find that 93% of the RF attributable to contrails is 895 due to the presence of natural clouds and that it decreases by 3% when accounting for contrail-contrail overlap. However, the magnitude of this effect is dependent on the optical thickness of the contrails, which remains highly uncertain (global estimations of average contrail optical depth can vary from ~0.065 to ~0.3).

Acknowledgments 900
The GEOS-FP meteorological data used in this study have been provided by the Global Modeling and Assimilation Office (GMAO) at the NASA Goddard Space Flight Center. This work was also partially funded by two separate grants: from the "la Caixa" Foundation through their Full Graduate Fellowship, and through NASA grant NNX14AT22A.

Code and Data availability
The codes and data used to produce this work are available from the authors upon request.

Author contributions
SB, RS, and SDE were responsible for the conceptual design of the study. FA, RS, and SDE provided 910 continuous review of the study progress and direction during execution. ISM performed all model simulations and wrote the manuscript. All authors contributed to manuscript editing. Qualitatively, the models behave similarly. Variation in longwave radiative forcing in response to changing optical properties is negligible in all three models, but our model consistently estimates a lower 1135 RFLW than is estimated by CoCiP. This error is maximized at low optical depths, reaching ~10%. The estimate from FL96 varies, agreeing more closely with CoCiP at low optical depths and more closely with our model at high optical depths.
Shortwave radiative forcing varies significantly with changes in optical properties in all three models. 1140 The range of asymmetry parameters simulated by our model results in a greater overall variation than is observed in the range of properties tested for FL96 or CoCiP. Qualitatively, the behavior of our model as the solar zenith angle (θ) increases matches that of CoCiP closely. RFSW increases slowly with θ, before reaching a peak between θ = 75° (high optical depths) and 88° (low optical depths). At values of θ beyond this peak, RFSW falls rapidly to 0. In FL96, the shape of the relationship is similar at all optical depths, 1145 with the minimum value occurring approximately at 75°. We also evaluate the difference in average value for each of the models, within the ranges of comparable microphysical properties. We obtain an average difference of less than 10% between CoCiP and our model, with the greatest error being ~20% at θ > 80°, expected from the here used two-stream approximation and mentioned in the limitations section. We find a greater average difference of ~40% between our model and FL96. 1150 To perform more quantitative analysis and comparison, we select a specific value of the relevant optical parameter for each model. For this purpose we choose an effective radius for ice of 45 µm. This is consistent with a natural ice cloud at an altitude of ~11 km, based on published parameterizations (Heymsfield and Platt, 1984;Corti and Peter, 2009;Lolli et al, 2017;Heymsfield, 2014). A prior analysis 1155 by Corti and Peter (2009) found that an asymmetry parameter g of 0.87 gave results which most closely matched those from FL96, and as such we use that value here. Key et al. (2002) also confirmed that this is consistent with solid columns of the mentioned size. Using the approach outlined earlier, this crystal size is represented in CoCiP using an effective radius of 45 µm and in FL96 using a generalized diameter of 46 µm. This specific single-contrail experiment results in differences in RFSW between our model and 1160 FL96 which are below 15%.

A2 Comparison of albedo effect on single contrail RFSW 1165
To evaluate the accuracy of our single contrail radiative forcing model, we simulate the effect of changes in surface albedo on single contrail shortwave radiative forcing and compare the results to both FL96 and CoCiP. Figure A2 shows the variation of RFSW with albedo in each model at three different optical depths.
As previously explained, CoCiP uses an effective radius of 45 μm and FL96 uses a generalized diameter of 46 μm, while our model is using an asymmetry parameter of g = 0.87. 1170 We observe the same qualitative behavior in all 3 models. Neglecting the already mentioned differences at high solar zenith angles (θ > 80°), FL96 and CoCiP quantitatively agree best with our model at low albedos (α < 0.5), with overall differences below 30%. Our model predicts a higher cooling impact at low solar zenith angles, and a lower cooling at high solar zenith angles. Maximum differences are found at 1175 high albedos ( α > 0.6) and high solar zenith angles (θ > 50°), where our model significantly underestimates RFSW, with differences of approximately 50%. However, these differences are less than 2 W/m 2 in absolute terms. The best agreement is found at an albedo of 0.3, the global Earth average albedo, with less than 10% difference. Finally, there are significant differences with CoCiP at low solar zenith angles and high albedos due to the forced negative sign of RFSW with that model. All percentage 1180 differences between the models are insensitive to changes in optical depth.

1185
To highlight these differences, we analyze one specific metric in all three models. Figure A3 shows the ratio of RFSW for an albedo of 0.5 to RFSW for an albedo of 0.3. In FL96, the ratio increases approximately linearly from ~0.1 to ~0.9 as the solar zenith angle increases from 0 to 90°. We also find a small increase in the ratio as a function of optical depth, although this falls with increasing solar zenith angle. 1190 cloud overlap with an effective decorrelation length (Barker, 2008). The simple expression of reflectance from Coakley and Chylek (1975), used in Corti and Peter model, allows us to develop our own formulation.

1215
In this section we develop the formulation for calculating shortwave radiative forcing for a 2-and 3-layers overlap and deduce a formulation applicable to an N-layers overlap. We start by recalling single contrail RFSW equation (Section B1.1), defined in main paper. We then develop the formulation for a 2-layers overlap (Section B1.2) and finish by extend the formulation to an N-layers overlap (Section B1.3), resulting in a simple formula easily applicable to our contrail coverage data. 1220

B1.1 Single layer RFSW
When evaluating shortwave radiative forcing of N infinite overlapping layers, we have to consider all the interactions between layers including reflectance and transmittance. We assume that cloud layers reflect 1225 shortwave radiation without diffusing it, whereas the Earth's surface diffuses incoming radiation in every direction (Corti and Peter 2009). Using these assumptions, we can decompose mathematically all the radiation interactions between layers.
As given in Section 2.1.1, the shortwave radiative forcing of a single contrail can be expressed as 1230 where S is the solar constant, is the Earth's surface albedo, t is the mean atmospheric transmittance, and R and R' are the direct and diffuse reflectances of the contrail. This expression can be rewritten as with (T, T') being the direct and diffuse transmittances (T = 1 -R, T' = 1 -R'). We can divide the shortwave to the contrail "albedo" for direct radiation, and is in this case simply R. The second, i2 can then be 1235 thought of as the contrail "albedo" for diffuse radiationin this case, ′ 1− ′ .

B1.2 Two-layer RFSW
Now consider a situation with two overlapping cloud layers, whose optical properties are fully captured 1240 by their individual reflectances ( 1 and 2 for direct, 1 ′ and 2 ′ for diffuse). Under the assumptions listed above, and ignoring edge effects, Fig. A1 diagrams the how one incoming unit of radiation (S = 1) will interact with these two layers.
1245 Figure B1. Decomposition of interactions between two cloud layers when receiving a single unit of shortwave radiation (S = 1). Sub-panels 1, 2, and 3 show three successive steps in the calculation as referred to in the text. R1: direct reflectance, T1: direct transmission (=1-R1); R'1: diffused reflectance, T'1: diffused transmission (1-R'1). Subscripts (1 and 2) indicate the layer number. = Earth albedo. Text color indicates the layer which most recently interacted with the radiation, with radiation from layer 1 in dark blue, from layer 2 in light blue, and from the 1250 Earth (reflected) in black. Subpanel 1 shows the initial interaction between incoming (direct) radiation and the upper contrail (contrail number 2). This contrail reflects a proportion R2 of the incoming direct light and transmits (allows to pass through) a fraction T2, equal to 1-R2. This would show all direct radiation interactions if 1255 there were no lower contrail, as light reflecting off the Earth's surface is assumed to be diffuse.
Subpanel 2 shows the full set of interactions for the incoming, direct, radiation when including both contrails. The fraction which passed through the upper contrail, T2, now undergoes an infinite number of reflections between contrails 1 and 2. On each reflection, some fraction (T1 and T2, respectively) of the 1260 reflected light passes through. This results in a geometric series, which can be summed to yield the total radiation which passes through the upper contrail (back to space) or lower contrail (towards the ground).
Ignoring these reflections, the radiation passing to the ground would be simply T1T2; the reflections increase this by a factor of 1 1− 11 2 , such that 1 2 1− 11 2 is the total radiation heading towards the surface.
The total which leaves upwards, back to space, is then 2 + α 11 2 2 1− 11 2 . 1265 Subpanels 3 then shows how diffuse radiation, reflecting off the Earth's surface, interacts with the system. As shown in Subpanel 2, the total direct radiation which reaches the ground is ( T 1 T 2 1− 11 2 ), of which only a fraction α is reflected back upwards as diffuse radiation. There are now two sets of infinite reflections to consider. The first is between the Earth and lower contrail, resulting in a geometric series which can 1270 be summed to 1 1−αR 1 ′ -now using the diffuse reflectance 1 ′ instead of the direct reflectance R1. The second is between the two contrails, and can be expressed using the effective "albedo" of the lower contrail 1 ′ (= 11 ′ + 12 ′ = 1 ′ + 2 ′ ). This geometric series can then be expressed as 1 1−α 1 ′ 2 ′ . From these equations, it becomes clear that the effect of additional contrails is to have additional "albedos", each of which modifies the total radiation which is either reflected to space or eventually absorbed by the Earth's surface 1275 (through repeated reflections).
The combination of the direct and diffuse radiation fluxes can then be seen in Subpanel 4; each upwards arrow from contrail 2 represents a separate component which will escape back to space. Adding these together and subtracting from the radiation which would be reflected to space under a clear-sky scenario 1280 (i.e. the Earth's albedo), the total shortwave RF can be summarized as where we have now combined the terms into two effective "albedos". These terms allow us to treat the combined layer pair as a if it were a single contrail. Specifically, we have 21 (= 2 + 11 2 2 1− 1 2 ) being the "albedo" of the layer pair to direct radiation, and 22 (= 12 2 2 ′ 1− 1 2 1 1− 1 ′ 2 ′ ) being the "albedo" of the layer pair to diffuse radiation. 1285

B1.3 N-layer RFSW
This approach extends from 2 to N layers by following the same mathematical logic (see Table B1), using as "albedo" values ( ) the direct and diffuse "albedos" of the (N-1) layers below the top one. 1290  This calculation means that we can collapse the effect of N different layers on shortwave radiation into the effect of a single, combined layer, as long as we know the direct and diffuse reflectances of each layer.
If we assume that all layers have the same optical properties (identical asymmetry parameter g and 1300 therefore identical optical parameter γ) we can simplify this further. Using the definition of R from the main text, we find that the direct albedos for 2 and 3 layers can be written as 21 = ( 1 + 2 )/ +( 1 + 2 )/ and 31 = ( 1 + 2 + 3 )/ +( 1 + 2 + 3 )/ . Extrapolating to an arbitrary N layers, we find that Therefore, the direct "albedo" from an N-layer overlap of similar layers (same optical properties) is equal to the direct reflectance of a single layer with the same total (summed) optical depth. The same logic can 1305 be applied to the diffuse albedo.
If the overlap occurs between layers of different optical properties, the same method can be applied as long as a single "effective" asymmetry parameter ge can be used for all layers. A method to find this parameter is derived below (Section B2). Once this parameter is known, RFSW for multiple overlapping 1310 contrails can be reduced to that for a single layer, i.e. B2 Derivation of weighted asymmetry parameter 1315 As outlined above, the calculation of shortwave radiative forcing for an N-layer overlap can be simplified significantly if a single, "effective" asymmetry parameter can be identified which characterizes the entire system. To calculate this effective optical parameter, we first determine what would be the effective optical parameter so that the direct radiation "albedos" are equal in both cases. We then show that matching this albedo is sufficient to ensure that the overall radiative forcing (accounting for both diffuse 1320 and direct radiation) matches between the "simplified" case and one in which each layer is treated independently.
In an N-layer overlap, a proportion N1 of incoming direct radiation is reflected. The single effective layer reflects radiation through the factor Re. The equality that must hold is then 1325 N1 = + (N−1)1 1 − (N−1)1 = (B9) As for eq. (A6), developing the expression of direct radiation albedo for a 2 and 3-layers overlap we obtain 21 = 1 2 / + 2 1 / 1 2 + 1 2 / + 2 1 / and 31 = 1 3 2 / + 2 3 1 / + 1 2 3 / 1 2 3 + 1 3 2 / + 2 3 1 / + 1 2 3 / . If we assume that an expression for the effective optical parameter of the entire layered system: If we use expression (B11) to calculate the effective diffuse radiation albedo ( ′ 1− ′ ) and expand each term, it results in the same formula for N2 as is shown in equation (A5). Since both the diffuse and direct albedos of the system are now matched, the total RFSW of the contrail layer system can be calculated by treating it as a single layer with the optical parameter shown in eq. (B11). 1335

B3 Variation of scattering with solar zenith angle
The objective of this section is to assess how the solar zenith angle (θ) affects the potential cooling impact from contrails. In Section 3.1.1 we stated that increasing the solar zenith angle θ also decreases the (net 1340 positive) contrail radiative forcing. This is because of an increase in shortwave cooling, since longwave radiation is not affected. Figure B2 shows how the total upscattered fraction of radiation is affected by changes in solar zenith angle. θ varies from 0 (noon) to 90° (sunset), moving anti-clockwise from the top left figure and shown as a black arrow. The dotted horizontal line represents the horizon. F and B represent downward (towards Earth) and upward (back to space) scattering. We assume that 90% of incident 1345 radiation is scattered forward, with 10% scattered backwards, representing the high forward scattering fraction of ice particles (high asymmetry parameter g). As the solar zenith angle increases, a greater fraction of the forward scattering peak is directed towards space (greater upscatter). This results in an increase cooling effect near sunrise or sunset compared to noon time.
1350 Figure C1 Components of response to a unit of incident light for 3 scenarios: no cloud/contrail, single layer, overlapping layers 1365 In order to explain this, we decompose the fraction of the incident SW radiation flux S reflected by layers of clouds or contrails into non-participating radiation ("NPR") and participating radiation ("PR") ( Fig.   D1). NPR is the light which is reflected into space from the upper contrail and therefore does not participate in scattering. In turn, it increases with Rc, which rises with optical depth. PR is the remaining outgoing shortwave radiation. Since all light included in PR was reflected or diffused, i.e. it has 1370 "participated" in scattering between the layer(s) and the Earth's surface, PR is driven by both direct reflectance Rc and diffused reflectance ′ . PR decreases with increasing optical depth. Finally, NR is the natural reflectance of light by the surface of the Earth, proportional to its albedo. We note that in the "clear sky" scenario, the total outgoing shortwave radiation is NR = ⍺S. With this decomposition, we can compute the shortwave RF of a single layer per unit of incoming radiation by comparing the outgoing shortwave radiation with no cloud (⍺) to that with a cloud layer https://doi.org/10.5194/acp-2020-181 Preprint. Discussion started: 5 May 2020 c Author(s) 2020. CC BY 4.0 License.