## Abstract

A new parameterized analytical model is presented to compute the instantaneous radiative forcing (RF) at the top of the atmosphere (TOA) produced by an additional thin contrail cirrus layer (called “contrail” below). The model calculates the RF using as input the outgoing longwave radiation and reflected solar radiation values at TOA for a contrail-free atmosphere, so that the model is applicable for both cloud-free and cloudy ambient atmospheres. Additional input includes the contrail temperature, contrail optical depth (at 550 nm), effective particle radius, particle habit, solar zenith angle, and the optical depth of cirrus above the contrail layer. The model parameters (5 for longwave and 10 for shortwave) are determined from least squares fits to calculations from the “libRadtran” radiative transfer model over a wide range of atmospheric and surface conditions. The correlation coefficient between model and calculations is larger than 98%. The analytical model is compared with published results, including a 1-yr simulation of global RF, and is found to agree well with previous studies. The fast analytical model is part of a larger modeling system to simulate contrail life cycles (“CoCiP”) and can allow for the rapid simulation of contrail cirrus RF over a wide range of meteorological conditions and for a given size-dependent habit mixture. Ambient clouds are shown to have large local impact on the net RF of contrails. Net RF of contrails may both increase and decrease and even change sign in the presence of higher-level cirrus, depending on solar zenith angle.

## 1. Introduction

Contrails are aircraft-induced cirrus clouds that may contribute to global warming (Minnis et al. 1999; Fahey et al. 1999; Burkhardt and Kärcher 2011). At high ambient humidity, contrails develop into contrail cirrus with properties similar to thin natural cirrus (Schumann 2002). Such contrails warm the Earth system during night, but may cool during day, in particular for large solar zenith angles (SZA), small particles, and particle habits with strong sideward scattering (Meerkötter et al. 1999; Myhre and Stordal 2001). Contrails tend to cool over dark and cold surfaces and tend to warm over bright and warm surfaces (Meerkötter et al. 1999). Contrails are often optically thin with optical depth (at 550 nm) of about 0.01–0.5 (Voigt et al. 2011). Thin cirrus in general tends to warm while, during daytime, thick cirrus tend to cool. Contrails may form above or below other clouds and are often associated with, or embedded in, thin cirrus (Sassen 1997; Immler et al. 2008). Ice particles in young contrails are typically small (1–20 *μ*m), and smaller than typical cirrus particles (Sassen 1997; Schumann 2002). Ice particles in contrails are nonspherical. The shape or habit of cirrus crystals varies in a not-well-known manner (Freudenthaler et al. 1996; Yang et al. 2010). Ice particles in young contrails are frozen droplets that may be described as droxtals (Roth and Frohn 1998). Droxtals differ from spheres in their scattering phase functions (Yang et al. 2003). With increasing age, the shape and size of contrail particles may approach that of natural cirrus. All these properties make simple radiation estimates difficult.

Several studies computed the changes in net downward irradiance (flux) or radiative forcing (RF) caused by additional thin cirrus and contrails (Stephens and Webster 1981; Meerkötter et al. 1999; Chen et al. 2000; Myhre et al. 2009; Rap et al. 2010; Frömming et al. 2011; Markowicz and Witek 2011). Rather complex models are required to accurately compute the scattering and absorption properties of ice particles (Key et al. 2002; Yang et al. 2003) and to compute the RF for given contrail properties, insolation, and atmospheric and surface properties (Fu and Liou 1993). The shortwave (SW) and longwave (LW) instantaneous RF by a contrail within a horizontally homogeneous atmosphere is the difference of the net downward flux values at the top of the atmosphere (TOA) with and without the contrail. The instantaneous RF is an approximation of the stratospherically adjusted RF commonly used in climate assessments (Hansen et al. 1997; Solomon et al. 2007). Contrails induce a small radiative cooling in the stratosphere (Meerkötter et al. 1999), causing only small differences between instantaneous and adjusted RF for contrails (Myhre and Stordal 2001; Marquart 2003; Hansen et al. 2005; Ponater et al. 2006). Although individual contrails differ from plane-parallel cloud layers (Schultz 1998; Gounou and Hogan 2007), contrails contribute most to RF when being long lasting, in which case they become rather wide. Therefore aged contrails are often approximated as thin plane layers in a horizontally homogeneous atmosphere (Minnis et al. 1999; Frömming et al. 2011). For global simulations with high temporal and spatial resolution and thousands of contrails, this task requires considerable amounts of computing time. Simplified relationships to estimate the RF from a layer of cirrus or aerosols have been suggested, but mainly for such layers in otherwise cloud-free atmospheres (Sobolev and Irvine 1975; Coakley and Chylek 1975; Charlson et al. 1992; Fu and Liou 1993; Baker 1997; Meerkötter et al. 1999; Kokhanovsky et al. 2005; Corti and Peter 2009).

Traditionally, RF of contrails and cirrus has been computed for given particle habit as a function of ice water path (Fu and Liou 1993; Yang et al. 2010), ice crystal size (Zhang et al. 1999), and size distribution (Mitchell et al. 2011). For given habit and particle size distribution, RF depends mainly on the optical depth *τ* (Ackerman et al. 1988). The optical depth *τ* is obtained by integrating the product of mass extinction coefficient *β* and ice (or liquid) water content over altitude (Hansen and Travis 1974). As shown in Key et al. (2002), for a wide set of particle size distributions and habits, *β* can be approximated as a function of the effective radius

where *V* is the particle volume and *A* is the mean projected particle cross-section area of the ensemble of particles in the cloud (Foot 1988; McFarquhar and Heymsfield 1998; Schumann et al. 2011). For constant *τ*, RF depends also, but less strongly, on the effective radius (Meerkötter et al. 1999).

In this paper, we develop a fast approximate model to estimate the LW and SW instantaneous RF at TOA for an added plane-parallel cirrus or contrail cloud. The method is designed to be part of an efficient tool for contrail cirrus prediction known as “CoCiP” (Schumann 2009, 2012). The tool simulates contrails resulting from a fleet of cruising aircraft, flight by flight, regionally or globally, for millions of contrail segments. CoCiP computes the optical depth and volume mean radius *r* of contrails for given ambient atmospheric conditions (Schumann 2012) and contains approximations to compute *r*_{eff} for given *r* (Schumann et al. 2011). The atmospheric conditions are taken from a host meteorological model (numerical weather prediction or global circulation climate model). Here, we consider the RF as a small disturbance relative to fluxes at TOA without contrails. The TOA fluxes are assumed to be known from the host model. These fluxes represent the effects of both cloudy and noncloudy ambient atmospheres on the radiation field.

The RF is computed for a set of ice particle habits (with random orientation) as a function of *τ* and *r*_{eff}. The RF results can be weighted according to empirically defined habit mixtures (Baum et al. 2005a; Schumann et al. 2011). The spherical approximation, which can be accurately calculated by Mie theory, is used here only for discussions and comparisons to previous results.

The contrail RF is computed using an approximate analytical model, basically a few equations, with some free parameters. The model equations (see section 3) are set up to approximate the functional dependence of the RF on variables characterizing the contrail and the ambient atmosphere. These functional dependencies are derived partly from previously suggested simplified relationships but mainly by selecting functions with as few free parameters as possible giving best approximations to forward calculations. The parameter values are determined by least squares fits to accurate forward calculations for a large set of test cases. The method design aims at limiting the model errors to values smaller than those from uncertain input values, in particular the ice particle habits and size distributions (Yang et al. 2010). The forward computations are performed using the libRadtran radiative transfer package by Mayer and Kylling (2005) (see section 2). The fit method and the resultant parameter values are reported in section 4. The model equations and the results are discussed, and validated for various applications, in clear and cloudy atmospheres, by comparison with a priori knowledge and published results in section 5.

## 2. Forward model data

Forward simulations have been performed for a large set of horizontally homogeneous atmospheric and surface conditions with “libRadtran” (Mayer and Kylling 2005). This code has been successfully validated in several model intercomparison studies and by comparison with observations (Van Weele et al. 2000; Mayer et al. 1997). The code libRadtran offers a flexible interface to set up the atmospheric and surface conditions as well as a choice of different radiative transfer equation solvers. The approach is as described previously (Krebs et al. 2007). For the present application, the irradiances are computed using the discrete ordinate solver by Stamnes et al. (1998), version 2.0, with six streams, which allows accurate simulations of irradiances. SW was calculated with the correlated-*k* distribution by Kato et al. (1999), which covers the wavelength range 0.24–4.6 *μ*m, and LW was calculated with the correlated-*k* distribution by Fu and Liou (1992), which starts from 4.54 *μ*m. Ice cloud single-scattering properties were parameterized using the double Henyey–Greenstein approximation (Key et al. 2002), except for droplets where Mie theory was used. To determine RF, the forward model is applied with and without an additional contrail layer.

To cover the variability of contrail RF with respect to contrail and ambient conditions, the dataset includes 4572 different atmosphere–surface cases for each of eight habit types. As sketched in Fig. 1, the cases represent different temperature and humidity profiles, over land and ocean, with and without upper-troposphere ice clouds and lower-level water clouds. Profiles of pressure, temperature, water vapor, ozone concentration, and other trace gases were taken from the Thermodynamic Initial Guess Retrieval (TIGR)-3 dataset (Chevallier et al. 1998). Although the surface emissivity in the thermal infrared is set to one, its variability is accounted for by the variability of the surface temperature. Reflection of radiation at land and ocean surfaces is represented by nonisotropic wavelength-dependent bidirectional reflectance distribution functions (BRDF). The ocean BRDF is parameterized as a function of the 10-m surface wind speed (Cox and Munk 1954; Nakajima and Tanaka 1983). For the land surface BRDF we used the Polarization and Directionality of the Earth’s Reflectances (POLDER)-1 BRDF database, issue 2.00 (Lacaze et al. 2002). The data are grouped into 17 International Geosphere Biosphere Project (IGBP) classes and 12 normalized difference vegetation index (NDVI) classes. The BRDFs were fitted using a three-parameter formula (Rahman et al. 1993) for randomly selected pixels from eight different Moderate Resolution Imaging Spectroradiometer (MODIS) scenes over woodland, grassland, snow, and desert. Each pixel has assigned a spectral albedo (determined by MODIS) as well as an IGBP class and an NDVI. The procedure covers a wide range of natural land surfaces.

For each case, 22 properties are varied randomly in certain ranges; 50% of all cases include a layer of water clouds, and 50% include a layer of ice clouds; 25% are cloud free, and 25% contain both water and ice clouds. The dataset applies for randomly oriented ensembles of particles, for seven different habit types (see Fig. 2): spheres (with 1–25-*μ*m radius), and solid columns, hollow columns, rough aggregates, rosettes, plates, and droxtals (with 10–45-*μ*m effective radius) (Baum et al. 2005a,b). The optical properties in the solar spectral range (volume extinction coefficient, single-scattering albedo, and asymmetry factor) of the particle ensembles are taken from Key et al. (2002), except for droxtals, which were provided by H. Gang and P. Yang (2009, personal communication). Optical properties for the thermal spectral range were taken from Yang et al. (2005). The data for droxtals and those in the thermal spectral range were processed exactly as the six habits described in Key et al. (2002), to obtain one consistent dataset. In addition, a so-called Myhre particle or flat particle with wavelength independent constant optical properties is used (Myhre et al. 2009; Markowicz and Witek 2011). These particles have an asymmetry factor of 0.8 for all wavelengths and a single-scattering albedo of 1.0 in the solar spectrum and 0.6 in the terrestrial spectrum. These values are typical for contrails (Strauss et al. 1997). Myhre particles are similar to solid hexagonal columns with fixed effective particle radius of about 32 *μ*m (Key et al. 2002).

In addition to changes in habit types, other properties were varied randomly in the dataset within the following ranges: cosine of SZA: 0.2–1.0 (SZA between 0° and 78°); surface temperature: atmospheric temperature at *z* = 0 km, ±10 K; ice cloud optical depth: 0–10; ice cloud bottom height: 6–10 km; ice cloud geometrical thickness: 0.5–2 km; water cloud optical depth: 5–50; water cloud effective radius: 5–15 *μ*m; water cloud bottom height: 1–2 km; water cloud geometrical thickness: 0.5–6 km; ocean surface wind speed (for ocean BRDF model): 1–15 m s^{−1}. Hence the dataset covers a wide range of atmospheric and surface conditions.

Contrails are included in these forward simulations with variable depth and altitude values, ice water contents, and particle habits. Contrail optical thickness is varied randomly from case to case over 0–2; contrail bottom height varies over 9–10 km; contrail geometrical thickness varies over 0.5–1 km. The atmospheric temperature at contrail midlevel varies between about 200 and 255 K. As a consequence, the optical depth of cirrus above contrails varies from 0 to 10.

Presently, the libRadtran code allows for only two cloud types per altitude bin, which we use to treat water and ice clouds. Therefore, the computations are performed assuming the particle habits and size distributions of contrail particles and ambient cirrus ice particles to be the same. This technical restriction, which may give reasonable results for given optical thickness, should be overcome in the future.

In total, the dataset contains 36 576 libRadtran cases with data for SW and LW RF. Figure 3 shows the computed RF versus various contrail or atmosphere properties. The panels show dependencies of the RF values on these properties, which form the basis for the model. In these data, the cooling in the SW dominates, mainly because of a large fraction of cases with large SZA, and nonspherical particle habits.

## 3. Model equations

The parameterization (the “model”) computes the changes in net downward LW and SW fluxes at TOA induced by an additional contrail layer (for 100% cover) for given Earth–atmosphere and contrail properties: In particular, the model uses the following input (W m^{−2}) as provided by meteorological models: outgoing longwave radiation (OLR) at TOA; reflected solar radiation (RSR) at TOA; solar direct radiation (SDR); *S*_{0}, the solar constant, with correction for the variable sun–Earth distance. Here, *μ* = cos(*θ*) = SDR/*S*_{0} defines the cosine of the solar zenith angle *θ* (SZA).

As can be seen from Fig. 3, RF_{LW} depends mainly on five independent properties, OLR, *T*: atmospheric temperature at the altitude of the contrail midlayer (K), *τ*: contrail optical depth at 550 nm, *τ _{c}*: optical depth at 550 nm of the cirrus above the contrails, and

*r*

_{eff}: the effective contrail ice particle radius (

*μ*m). All of them are incorporated into the model.

For each habit class, we use the following empirical equation to compute RF_{LW}:

The five (positive) model parameters for LW are *T*_{0}, *k _{T}*,

*δ*,

_{τ}*δ*

_{lc}, and

*δ*

_{lr}, as explained below. The number of parameters is the minimum needed for five independent properties. The parameter values, as determined in section 4, are listed in Table 1.

The first term in Eq. (2) estimates the change in outgoing LW radiation flux assuming opaque contrails. The LW flux from space is small and hence neglected in this approximation. The contrail temperature *T* is usually lower than the brightness temperature (related to OLR by the Stefan–Boltzmann law) of the atmosphere without contrails. The linear approximation with free parameters *K _{T}* and

*T*

_{0}has been found to approximate this flux change better than any variant of the Stefan–Boltzmann law. The second factor accounts for the effective emissivity of the contrail: Here the absorption optical thickness in the thermal IR is the relevant quantity, which is typically only 50% of the optical thickness at 550 nm. This is considered by the product of a parameter

*δ*and a factor

_{τ}*F*

_{LW}accounting for the size dependence of the optical properties (extinction relative to scattering and absorption):

The last factor approximates the reduction of the OLR at the contrail level due to a cirrus with optical depth *τ _{c}* above the contrail:

The LW RF is constrained to positive values, so that contrails at low altitudes with high ambient temperatures have zero contribution to RF_{LW}.

From Fig. 3 we see that RF_{SW}, for given solar direct irradiance SDR, depends on five independent properties: *τ*, *μ* = cos*θ*, effective albedo of the Earth–atmosphere system *A*_{eff} = RSR/SDR, *r*_{eff}, and *τ _{c}*. This dependence is parameterized by

The 10 model parameters for SW are *t _{A}*, Γ,

*γ*,

*A*,

_{μ}*B*,

_{μ}*C*,

_{μ}*F*,

_{r}*δ*

_{sr},

*δ*

_{sc}, and . In view of the number of independent properties (5) and the complexity of light scattering, the number of free parameters is still small. The parameters are defined to be positive, with values as determined in section 4 and listed in Table 1.

The first factor in Eq. (5) is the incoming solar radiation. The second factor describes the dependence on the effective albedo *A*_{eff}. The third factor can be interpreted as the albedo of the contrail:

Its values depends on effective optical depth values:

which are functions of SZA and particle size,

reflectances

and an empirical function accounting for SZA-dependent sideward scattering

Finally,

accounts for the optical depth *τ _{c}* of the cirrus above the contrail and its effective value

*τ*

_{c}_{,eff}=

*τ*/

_{c}*μ*. The model computes negative RF

_{SW}values; small positive values, which may occur rarely in absorbing atmospheres (Meerkötter et al. 1999), are not predicted by this model.

The model computes the LW and SW contributions to RF as a sum of the contributions from various habit classes:

Here *G*_{habit} = *G*_{habit}(*r*) is the contribution of each habit to RF. For contrails with effective particle radius below 30 *μ*m, we assume a size-dependent habit mixture consisting of droxtals, solid columns, and 3D bullet rosettes (Schumann et al. 2011). Larger contrail particles are assumed to have the same habits as natural cirrus (Baum et al. 2005a). The habit model as a function of volume mean radius *r* is given in Table 2 of Schumann et al. (2011). For each habit, *r*_{eff} is computed for given *r* with relationships as given in Table 1 of that reference.

The model has been coded in FORTRAN 77. A single evaluation of both the SW and LW contributions together requires about 22 *μ*s on a modern laptop, about a factor 10^{5} less than in libRadtran. This makes extensive contrail simulations feasible.

## 4. Model parameters

The model parameters (see Table 1) have been computed by a least square fit where the sum

over all data *d _{i}* and model results

*f*

_{i}=

*f*

_{i}(

**x**) (

*i*= 1, 2, … ,

*N*) assumes its minimum for variations in the free parameters that compose the vector

**x**. Here,

*S*

^{1/2}is the root-mean-square (rms) error. The weights

*G*> 0 may account for the frequency distribution of independent properties. Here, we set

_{i}*G*= 1 + 1/(

_{i}*τ*+ 0.1) to give cases with low contrail optical depth

_{i}*τ*a higher importance in the fits.

_{i}Several algorithms exist to find the minimum of *S*(**x**); see, for example, Yang et al. (2000). In principle, the minimum is characterized by zero derivatives of *S* with respect to **x**. In a certain vicinity around the optimum at **x**, the function *S*(**x**) exhibits several local minima, some of which with equal value, so that a unique absolute minimum and well-defined derivatives may not exist. The diameter of this vicinity defines an uncertainty range for the parameters.

Here we use a simple but reliable, robust, and sufficiently efficient search method that determines the minimum without explicit knowledge of the derivatives of *S*, overcomes local minima, and determines the uncertainty range of the results. The search method starts from an arbitrary initial guess **x**_{0}, with *S*(**x**_{0}) > 0, and initial steps with sizes on the order of the magnitude of **x**_{0}. The method searches for lower approximation errors in the neighborhood of the first guess. For this purpose, the method loops over all components *i* of the vector **x** and computes for . The search first considers positive changes, *s* = 1, and, if no lower *S* value is found there, negative ones, *s* = −1. If , we replace **x**_{0} with and start the search again. Otherwise, the search for smaller *S* values continues for the next component *i*. The search for given Δ**x** ends when *S* is a minimum relative to all neighboring values. Thereafter the search is repeated with progressively smaller step sizes Δ**x** until the step size drops below a prescribed small value *ε* of order 10^{−6}. To overcome local minima, the search is repeated several times starting from the last solution, again first with large initial steps Δ**x**. At the end, a search over the vicinity of the solution **x** with the smallest step sizes identifies local minima in the neighborhood. The diameter of the domain containing such local minima defines the uncertainty range of the model parameters for the given database. Largest uncertainties (8%) are found for *δ*_{sr}, and smallest are found (5 × 10^{−4}%) for *T*_{0}.

Table 1 lists the fitted parameter values. Figure 4 shows the approximation errors in relation to the SW and LW and net (LW + SW) RF values. The errors are significant but are small relative to the various uncertainties introduced. Figure 5 shows that the model and the libRadtran results for SW and LW RF correlate at better than 99%, with mean biases below 2%. For RF_{LW} (see Fig. 5) the correlation coefficient is *r* = 0.995, and the rms error is about 2.1 W m^{−2}. For RF_{SW}, the correlation coefficient is 0.991 with rms error of 4.0 W m^{−2}. On average over all habits, the relative errors are 7.1% for LW and 10.6% for SW fluxes.

The fit results depend on the range of the selected properties in the libRadtran database. The robustness of the fit has been investigated by repeating the fit calculations with subsets of the test dataset. Table 2 shows the percentage rms errors of the model for contrails containing solid hexagonal columns. The first line lists the results of the present model with parameter values as given in Table 1. The model fits the LW forcing better than the SW forcing, but both with errors below 10%. The same model applied to only cloud-free test cases gives slightly smaller errors. A separate fit (“refitted”) especially for cloud-free cases reduces the errors only slightly. For cloudy atmospheres, the RF magnitude is smaller and hence the relative errors are larger. The model approximates both cloudy and cloud-free atmospheres with similar absolute accuracy. The model error increases slowly with *τ* (see Fig. 4). Tests for thick contrails (1 ≤ *τ* ≤ 2) showed slightly smaller relative errors than for the whole range of *τ* values. For a case restricted to thinner contrails (*τ* < 1), the absolute rms errors are also significantly smaller. Hence the model is suited both for optically thin and thick contrails.

The different model terms are all required to reach high accuracy, but to a different degree. The LW rms errors are most sensitive to the value of *T*_{0}. When the parameters *T*_{0}, *k _{T}*,

*δ*, and

_{T}*δ*

_{lc}are individually varied by 10%, and the other parameters refitted to minimum rms errors, the rms errors increase by 60.9%, 22.6%, 21.8%, 0.08%, and 0.005%, respectively. For the SW parameters, the largest changes (70.8%) occur for

*t*and the smallest (0.08%) for

_{A}*δ*

_{sc}. The latter sensitivity is small partly because cases with cirrus above contrails are rare in the dataset and partly because the RF in cloudy atmospheres is smaller compared to the cloud-free case anyway.

The remaining deviations between the model and the libRadtran data are due to variations of atmosphere–surface properties not explicitly covered by the model. In particular, the model does not explicitly account for contrail altitude and thickness, humidity, presence of clouds below the contrails, and variations in surface properties. The dependence of RF on geometrical contrail depth is negligible in the libRadtran results.

The largest relative errors occur for SW cloudy land cases (18%), and the largest absolute error (7.2 W m^{−2}) for SW cloud-free land, mainly because of the anisotropic BRDF of the Earth surface. Far smaller errors are found over ocean (relative errors for SW cloud-free ocean 2.5%). Hence, the errors could be reduced significantly by using separate sets of parameter values for cloudy and cloud-free atmospheres and for land and ocean surfaces, but this would complicate the model and require further input for applications.

## 5. Discussion

### a. Comparison for idealized cloud-free cases

We discuss the model in comparison with previous results, first for cloud-free atmospheres. Figure 6 compares the RF values for a contrail with midlevel at 10.5-km altitude containing spherical particles with a previous “benchmark” result for a cloud-free summer atmosphere above a Lambertian surface with a surface albedo of 0.2, obtained with the matrix operator model (Meerkötter et al. 1999). This benchmark has also been used elsewhere (Myhre and Stordal 2001; Stuber and Forster 2007; Rap et al. 2010; Markowicz and Witek 2011). As input for our model in this comparison, we need to know the OLR and RSR values. Both values, see legend of Fig. 6, have been provided from the previously published computations (R. Meerkötter 2009, personal communication). Here, RSR and SDR are functions of the SZA, for 45° latitude, 21 June.

We find excellent agreement for daily mean LW and SW RF (evaluated with 1-min time resolution). As is well known (Stephens and Webster 1981), the cloud emissivity and hence RF_{LW} saturates faster with increasing *τ* than cloud albedo and RF_{SW} (see Figs. 3 and 6). The model simulates an effective contrail albedo *α _{c}*(

*μ*,

*τ*), which increases nearly linearly with

*τ*for

*τ*< 3 and SZA < 80°. This albedo is smallest for spherical particles and largest for solid hexagonal particles (varying from 0.13 to 0.25 for

*τ*= 0.3, SZA = 0).

The coefficient Γ in the reflectance *R _{C}* of a nonabsorbing cloud layer is related to (1 −

*g*), where

*g*is the asymmetry factor, a measure for the fraction of light scattered forward by the cloud particles (Coakley and Chylek 1975; Baker 1997; Corti and Peter 2009). Consistent with this interpretation, Γ is smallest for largest

*g*, that is, for spheres compared to other habits (see Table 1).

The RF_{SW} tends to zero for *A*_{eff} → *t _{A}* < 1 (see Fig. 3). The quantity

*t*is the albedo of an atmosphere in which an additional cirrus cloud does not enhance the albedo. This value is less than 1 even over ideally reflecting surfaces because of absorption in the atmosphere. For Myhre particles, with zero SW absorption,

_{A}*t*is largest compared to other habits and close to 1. The quadratic dependence on (

_{A}*t*−

_{A}*A*

_{eff}) accounts for multiple reflections between the contrail and the surface (Charlson et al. 1992). Tests with a term

*t*(1 −

_{A}*A*

_{eff}) (Corti and Peter 2009) instead of (

*t*−

_{A}*A*

_{eff})

^{2}in Eq. (5) turned out to be less accurate. Tests with a variable exponent instead of the constant 2 showed that the value 2 is close to optimal.

For comparison, Table 2 shows the performance of the parameterization of Corti and Peter (2009) for cloud-free cases as a function of surface temperature and surface albedo, first with the published parameter values in comparison to the libRadtran results, and second with parameter values fitted to the present test data. In both versions, the present model shows significantly smaller errors.

The libRadtran data reveal a decrease of the effective albedo *A*_{eff} with SZA. This suggests an alternative model with a SZA-dependent *t _{A}* values and a further model parameter. However, tests have shown that this does not reduce the errors significantly. We also tested model variants in which we replaced terms with exponentials 1 − exp(−

*αx*) with

*αx*/(1 +

*αx*), and found that the given variant is the one with lowest rms errors. An exponential dependence on optical depth is expected for thin cirrus (Sobolev and Irvine 1975; Kokhanovsky et al. 2005).

The LW model parameterizes absorption and emission by the contrail for given contrail temperature *T* linearly by *k _{T}*(

*T*−

*T*

_{0}). The value of

*T*

_{0}is found close to 150 K, without obvious physical meaning. This relationship could be interpreted as a linearization of a dependence on . A fit with such a function, as suggested by Corti and Peter (2009), gives an exponent of about 3. This value is considerably smaller than Stefan–Boltzmann’s temperature exponent 4 for radiation from a blackbody in vacuum because of absorbers above the contrail (Corti and Peter 2009). Hence, the coefficient

*k*may depend on the water vapor column above contrails. In fact, as also shown by the libRadtran data, OLR versus surface temperature is best approximated with an even smaller exponent (about 2), because of absorbing clouds and water vapor in the atmosphere above the surface. However, somewhat to our surprise, the linear fit performs best for the present set of forward calculations.

_{T}The RF_{LW} value increases linearly for small values of *τ* but gets saturated quickly as modeled by the emissivity [1 − exp(−*δ _{τ}τ*)]. The value of

*δ*is habit dependent. It is smallest (0.68) for rough aggregates and largest for spheres, inversely proportional to the particle size effective for thermal absorption. For solid hexagonal columns, the value of

_{τ}*δ*(0.81) is close to the values 0.79 and 0.75 derived before by Fu and Liou (1993) and Corti and Peter (2009), respectively. This agreement is satisfactory in view of different radiation transfer models and atmospheres used for calibration in the studies cited.

_{τ}Figure 7 shows that our model approximates the SZA dependence as reported in Meerkötter et al. (1999) to reasonable accuracy. Part of the differences are caused by the surface reflectance function, which was Lambertian in the benchmark computations but include more realistic anisotropic cases in our fit. For SW forcing, the SZA dependence is very important. The effective optical depth *τ*_{eff} = *τ*/*μ* causes a monotonic decrease of SW forcing for growing SZA. It does not explain the maximum in RF_{SW} around SZA 70°–80° as shown in Fig. 8. Here, inclusion of the function *F _{μ}*(

*μ*) was essential to describe this strong dependence of RF

_{SW}on SZA. This dependence is known to be stronger for thin cirrus than for thick ones (Coakley and Chylek 1975), and this fact is accounted for by multiplying

*F*with the reflectance , which is largest for small

_{μ}*τ*and small SZA.

The RF by contrails depends strongly on the ice particle habit, as can be seen from Fig. 3, and specifically from Figs. 8 and 9. The latter results have been computed for conditions similar to those considered in Markowicz and Witek (2011). Here we had to estimate their OLR and RSR values (see Fig. 3). Nevertheless, the range of values and the shape of the curves are very similar to those shown in that paper.

Figure 8 stresses the importance of particle habits for daytime contrail RF. Whereas spheres may cause a net warming even for quite thick contrails, the other habits may cause a net cooling of the atmosphere during daytime even for rather thin contrails because of far stronger sideward scattering causing (up to 100%) stronger RF_{SW}. Also RF_{LW} depends on particle habit. The LW difference between spheres and rough aggregates amounts to 30%. For other habits the differences remain below 20%.

RF depends also on the effective radius of the contrail particles (see Fig. 3). Figure 10 shows the size dependence of RF for constant *τ* and constant SZA for the range of sizes covered by the data. This dependence reflects the size dependence of the ratio of single-scattering properties of the particles relative to the extinction determining the optical depth (Key et al. 2002; Yang et al. 2005). The radius dependence is stronger for SW than for LW radiation (Baum et al. 2005b). The increase with radius approximated by the functions *F*_{LW,SW}(*r*_{eff}) is caused by increasing absorption and forward scattering. The absorption coefficient of ice at 10-*μ*m wavelength is about 0.06 *μ*m^{−1} corresponding to an absorption length of 16 *μ*m (Warren 1984). Therefore, particles with radius in the order of 10 *μ*m or more are basically opaque to radiation. For smaller radius, RF_{LW} is smaller. Droxtals behave similar to spheres in the LW range (because of nearly equal *r*/*r*_{eff} and volume/area ratios) but show stronger RF for SW because of more sideward scattering caused by the faceted surface. The size dependence of the RF values for plates is different from those of other habits because of smallest reflectance and absorption efficiency (Key et al. 2002; Yang et al. 2005) for the same effective radius as a consequence of their special geometry with rather small *r*_{eff}/*r* ratio (Schumann et al. 2011). For the Myhre particles, the size dependence vanishes because of constant asymmetry factor, zero absorption in SW, and constant absorption in LW.

For constant *τ* and constant SZA, Fig. 10 shows that the net RF is largest for large radii *r*_{eff}. However, for constant ice water path (IWP), the optical depth increases with decreasing effective radius (at least for particle diameters large compared to the solar wavelength), so that the net RF is largest for small particles. This is consistent with previous studies (Hansen and Travis 1974; Fu and Liou 1993; Zhang et al. 1999).

### b. Comparisons for idealized contrails in a global cloudy atmosphere

As a more realistic test, we applied our model to the intercomparison case of Myhre et al. (2009). This test case (see also Rap et al. 2010; Frömming et al. 2011), compares instantaneous TOA RF from various global models with state-of-the art radiation transfer models and realistic distributions of temperature, clouds, and surface properties, for one year of meteorological data. In the intercomparison, each modeling group used their own meteorological data (from a numerical weather prediction or a global circulation climate model). Annual mean values of the forcing results are used for comparisons. The test considers a 1% homogeneous contrail cover with fixed contrail properties around 10.5-km altitude and optical depth of 0.3. To ensure that differences in single-scattering properties of contrail particles do not influence the intercomparison, Myhre particles with wavelength-independent properties were specified, see section 2. The comparisons were performed both for clear sky and an “all sky” case. Here we compare to the results for the all-sky case, that is, for a uniform contrail layer in a cloudy atmosphere.

The present model is applied using forecast data (3 hourly) for pressure, temperature, cirrus ice water content IWC_{C}, and radiative fluxes (OLR, RSR, SDR) at TOA for the year 2006 with 1° horizontal resolution from the Integrated Forecast System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF; see http://www.ecmwf.int/research/ifsdocs/). The optical depth of cirrus above contrails is computed from these data for given cirrus-extinction *β _{C}* = 3

*Q*

_{Cext}IWC

_{C}/(4

*ρ*

_{ice}

*r*

_{effC}), ice bulk density

*ρ*

_{ice}= 917 kg m

^{−3}, and cirrus particle extinction efficiency

*Q*

_{Cext}= 2. The effective radius is determined as a function of temperature and ice water content of the ECMWF cirrus (Sun 2001; Sun and Rikus 1999). In this application we assume that the OLR and RSR values account for broken cloud cover and assume that the contrails are randomly distributed.

Figure 11 shows the SW, LW, and net annual mean RF fields computed with the present model, averaging over 3-hourly results. The results from our model are similar to those of five models shown in Myhre et al. (2009). The geographical distribution exhibits, as expected, a maximum negative SW forcing over dark surfaces (cloud-free ocean areas), while the SW forcing is small over snow (Arctic and Antarctic) and high clouds (e.g., in the intertropical convergence zone). Maximum LW forcing is computed for a contrail layer over the warm and cloud-free deserts (in North Africa and Australia) and in the subsidence branches of the Hadley cell over the subtropical oceans. The net results differ more strongly between the models. Two out of the five models find slightly negative net values regionally. Our model gives positive mean values everywhere for this application, but the values are close to zero over the maritime continent in the tropical western Pacific, the Himalaya, and parts of the Antarctic.

The mean values of this model and of the other models [from Fig. 5 in Myhre et al. (2009) for all sky] are listed in Table 3. The model acronyms are explained in the reference paper. The statistical uncertainty of the annual mean values obtained from daily mean values is 2%. We see that the global mean RF values computed with the present model for this intercomparison case agree very well with the results from the other models. It deviates by less than 10% from the results of the first model listed in this table, which uses meteorology data and radiative transfer models similar to those used here.

### c. Contrails above, within, or below clouds

The RF by contrails changes considerably when forming above, within, or below other clouds. Low-level clouds generally cause high RSR (high albedo) and low OLR. Hence they reduce both the negative SW forcing and the positive LW forcing. Figure 3f shows that the changes in SW RF dominate, so that low-level water clouds generally cause an increase in net RF of contrails relative to cloud-free cases. Figure 3e shows that contrails contribute considerably to RF even when forming in or below ambient cirrus of optical depth up to about 10.

The present model accounts for the change in RF by reduced RSR and OLR, which enter the model equations directly. The additional change of contrail RF by cirrus above the contrail is not very strong (Fig. 3e) relative to the clear-sky case. It is approximated as a function of cirrus optical depth *τ _{c}* in the factors

*E*

_{LW,SW}[Eqs. (4) and (11)]. Figure 12 shows the correlation between the libRadtran results and the model for cases with cirrus above the contrails. The correlation is not quite as good as for the whole dataset, but still reasonable (correlation coefficients larger than 0.98 both for SW and LW). Because of the weak influence of

*τ*, the correlation coefficients depend only weakly on the model parameters in

_{c}*E*

_{LW,SW}. Note that the fit is limited to optical depth of cirrus above contrails between 0 and 10. This range should be sufficient for most applications. In the global intercomparison case (section 5b), the cirrus above the contrail is thin (

*τ*< 3) because of low ice water content due to low ambient temperatures above contrails at 10.5-km altitude. Hence, the model may be applied also for contrails above or below other clouds.

_{c}The libRadtran results reveal that RF_{SW} gets stronger (more negative) for small SZA but weaker for large SZA. This effect is reflected by the different signs of the terms with *τ _{c}* and

*τ*

_{c}_{,eff}in

*E*

_{SW}, Eq. (11). For example, for solid column particles and

*τ*= 3, the factor values are

_{c}*E*

_{LW}= 0.75 and

*E*

_{SW}= 1.15 0.34 for SZA= 20° and 75°, respectively. Hence, high-level cirrus may both increase and decrease the SW RF depending on SZA. The SZA-dependent changes are caused by either larger (for small SZA) or smaller (for large SZA) downward SW radiances reaching the contrail because of sideward scattering of solar light by the cirrus particles. The LW RF of contrails below high-level cirrus is always smaller than below clear air. Hence, the local net RF of contrails may change sign when forming in the presence of cirrus instead of clear air.

## 6. Conclusions

An analytical parametric model has been presented that allows fast computations of the radiative forcing by a contrail or an added cirrus layer. Unlike earlier studies, the RF is computed not as a function of surface properties (albedo and temperature) but as a function of TOA reflected shortwave and outgoing longwave radiation. These fluxes are available from numerical weather prediction or climate model results, and may be measured from space. This makes the model applicable both to cloud-free and cloudy atmospheres. The model accounts for dependence of RF on contrail temperature, optical depth, effective particle radius, solar zenith angle, and the optical depth of cirrus above the contrail. In particular, the model captures the RF_{SW} minimum at large solar zenith angles. The model includes eight habit types. RF results may be weighted according to empirically defined habit mixtures, but data are missing to test this approach. The model approximates forward calculations with an accurate radiation transfer model with relative errors of about 7% for LW and 11% for SW RF. In view of other uncertainties—for example, differences in RF for different habits of the order 100% as exhibited in Figs. 8 and 9, and similarly large differences between various global model results for a well-defined intercomparison case (Table 3)—this appears to be acceptable for practical applications. The main advantage of the model is its simplicity, which makes the simulation of millions of contrail cases possible within seconds of computing time.

The model has been validated for various test cases. In agreement with other studies we find a strong sensitivity of contrail RF to solar zenith angle and particle habits. For constant optical depth, the RF is slightly dependent on the particle size because the ratio of extinction to scattering and absorption changes with particle size. In addition, we have discussed the impact of ambient cirrus on contrail RF. Low-level clouds below the contrail generally increase the local net contrail RF by stronger reduction of the SW RF relative to the LW part. Contrails contribute considerably to SW RF even when forming in or below ambient cirrus of optical depth up to about 10. The net RF of contrails may both increase and decrease and even change sign in the presence of cirrus depending on solar zenith angle.

The model may be improved in future studies. In particular, the habit mixture should be included already when setting up the optical properties of the particle ensembles in the forward model calculations. This would also further reduce the computing time. The forward model should be extended to include contrails with habit and size distributions different from those in ambient cirrus. Effective radius values below 10 *μ*m should be included for all habits occurring in thin cirrus and contrails, in particular droxtals. Finally, other functional dependencies may be found, possibly including other atmosphere–surface or contrail properties approximating the forward simulation results with smaller deviations. Atmospheres may be taken from the output of numerical weather prediction models for meteorological and air traffic conditions where contrails occur, which may then also allow to account for the impact of the water vapor column changes above contrails. The present study is restricted to particles with random orientation; other orientations may be of interest as well. Finally, if three-dimensional forward simulations would be possible for a large set of parameter values, one could also account for more realistic contrail geometries.

## Acknowledgments

The POLDER-1 BRDF database issue 2.00 has been produced by POSTEL from an original algorithm developed by Noveltis. The data are from CNES POLDER-1 on board NASDA *ADEOS-1*. We thank Drs. Luca Bugliaro, Ulrich Hamann, Alexander Kokhanovsky, Ralf Meerkötter, Gunnar Myhre, Michael Ponater, Margarita Vázquez-Navarro, Christiane Voigt, and Eleonora P. Zege for input and valuable suggestions. We also acknowledge the use of ECMWF data, with permission from the German Weather Service. This work contributes to a DLR project CATS and was partly funded by a German BMBF project UFO and an EU-project REACT4C.

## REFERENCES

## Footnotes

Additional affiliation: Meteorologisches Institut, Ludwig-Maximilians-Universität, Munich, Germany.