Generated using the official AMS LATEX template—two-column layout. FOR AUTHOR USE ONLY, NOT FOR SUBMISSION!
JOURNAL

OF THEATMOSPHERIC SCIENCES

Entrapment: An important mechanism
to explain the shortwave 3D radiative effect of clouds
ROBIN J. H OGAN∗ , M ARK D. F IELDING
European Centre for Medium-Range Weather Forecasts, Reading, UK

H OWARD W. BARKER
Environment and Climate Change Canada, Toronto, Canada

NAJDA V ILLEFRANQUE
CNRM, Météo-France/CNRS, Toulouse, France, and
LAPLACE, Université de Toulouse, Toulouse, France

S OPHIA A. K. S CH ÄFER
Meteorological Institute, Ludwig-Maximilians-Universität München, Munich, Germany

ABSTRACT
Several mechanisms have previously been proposed to explain differences between the shortwave reflectance of realistic cloud scenes computed using the 1D Independent Column Approximation (ICA) and3D solutions of the radiative transfer equation. When the sun is low in the sky, interception of sunlight bycloud sides tends to increase reflectance relative to ICA estimates that neglect this effect. When the sun ishigh, 3D radiative transfer tends to make clouds less reflective, which we argue is explained by the mechanism of “entrapment” whereby horizontal transport of radiation beneath a cloud layer increases the chances,
relative to the ICA, of light being absorbed by cloud or the surface. It is especially important for multi-layered
cloud scenes. We describe modifications to the previously described Speedy Algorithm for Radiative Transfer
through Cloud Sides (SPARTACUS) to represent different entrapment assumptions, and test their impact on
65 contrasting scenes from a cloud-resolving model. When entrapment is represented explicitly via a calculation of the mean horizontal distance traveled by reflected light, SPARTACUS predicts a mean “3D radiative
effect” (the difference in top-of-atmosphere irradiances between 3D and ICA calculations) of 8.1 W m−2
for overhead sun. This is within 2% of broadband Monte Carlo calculations on the same scenes. The importance of entrapment is highlighted by the finding that the extreme assumptions in SPARTACUS of “zeroentrapment” and “maximum entrapment” lead to corresponding mean 3D radiative effects of 1.7 W m−2 and19.6 W m−2 , respectively.

1. Introductionof the Speedy Algorithm For Radiative Transfer throughCloud Sides (SPARTACUS; Hogan et al. 2016; Schäfer
et al. 2016). This solver is now available as an option in
the radiation scheme of the forecast model used by theEuropean Centre for Medium-Range Weather Forecasts,
ECMWF (Hogan and Bozzo 2018).
In the shortwave, the main effect of transport throughcloud sides is “side illumination”: the enhanced interception of direct sunlight when the sun is low in the sky,which increases the reflectance of the scene and equivalently the magnitude of the cloud radiative effect (CRE).
However, Barker et al. (2015) reported Monte Carlo calculations showing that in realistic cloud scenes, the effect of
introducing 3D transport was more typically to reduce the

A key challenge in atmospheric modeling for bothweather and climate prediction is to improve the interaction of clouds with solar and thermal-infrared radiation. While the representation of sub-grid cloud horizontal structure and vertical overlap is often now quite
sophisticated (e.g. Pincus et al. 2003; Hill et al. 2015;
Di Giuseppe and Tompkins 2015), a process missing fromall operational models is the horizontal transport of radiation within gridboxes. This was characterized by Hogan
and Shonk (2013) as entirely associated with flow of radiation through cloud sides, and led to the development∗ Corresponding author address: Robin J. Hogan, ECMWF, Shinfield Park, Reading, RG2 9AX, UK.
E-mail: r.j.hogan@ecmwf.int

Generated using v4.3.2 of the AMS LATEX template

1

2

JOURNAL

OF

THE

ATMOSPHERIC SCIENCES

magnitude of the CRE, particularly when the sun is highin the sky.
Várnai and Davies (1999) characterized 3D solar radiative effects in single-layer cloud scenes in terms of fourmechanisms, two of which reduce the reflectance of acloudy scene and are therefore candidates to explain thisbehavior. Their “downward escape” mechanism explains
how forward-scattered sunlight inside a cloud has a chanceto escape through the side of a cloud and reach the surface, whereas in the corresponding ICA case it would remain within the cloud and have more chance of being scattered back to space. Welch and Wielicki (1984), Hogan
and Shonk (2013) and Barker et al. (2016) argued that thisescape mechanism explains why 3D effects reduce the reflectance of cumulus, stratocumulus, aircraft contrails andstochastically generated cloud fields, for high-sun conditions. This process is represented by SPARTACUS, but inthis paper we present evidence to show that it is not significant enough to explain the results of Barker et al. (2015),
which were for a wide range of realistic and often multilayered cloud scenes.
The second candidate mechanism from Várnai andDavies (1999) is “upward trapping”, which incorporates
all light rays that (i) are reflected back to space in the
ICA case but not when 3D transport is included, and (ii)have a longer path length in 3D than ICA. Their diagram
to explain how this mechanism typically acts in singlelayer cloud scenes depicted an upward traveling light ray
passing horizontally through the side of a cloud above, a
process that is already represented by SPARTACUS. However, it is also possible, particularly in multilayered scenes,
for trapping to occur without any transport through cloudsides, but rather as a consequence of horizontal transport
entirely within a clear or cloudy region and the upwardreflected ray then intercepting the base of a cloud above.
This process is not explicitly handled by SPARTACUS. It
was alluded to by Barker and Davies (1992) who considered idealized single-layer clouds over a reflective surface,
but not studied in detail for realistic multilayered cloudscenes.
In this paper, we seek to quantify the importance of thismechanism, which we refer to as “entrapment”. In section
2, we describe it in more detail and present a simple mathematical example to illustrate how it reduces the scene reflectance. In section 3, we describe how the limits of zeroand maximum entrapment may be represented in SPARTACUS. This is followed by two sections on the more
complex “explicit” entrapment calculations: section 4 describes how we estimate the horizontal distance traveledby reflected radiation, with validation against monochromatic Monte Carlo simulations, while section 5 describes
how the distance traveled is used to compute how much
entrapment occurs, accounting for the fractal nature ofclouds. Readers uninterested in the internal workings of

SPARTACUS may wish to skip sections 3–5. Then in section 6, estimates of the broadband shortwave 3D radiative effect by the new SPARTACUS solver are evaluatedby comparing to Monte Carlo calculations performed on
65 diverse high-resolution scenes from a cloud-resolving
model.
2. The concept of entrapment
The schematic in Fig. 1 illustrates how entrapment canchange the reflectance of a cloud scene. Panel a depicts the
behavior assumed in the ICA, in which horizontal transport is ignored: incoming solar radiation scattered upwardby the first cloud layer it encounters is likely to escape tospace since it passes back through the same clear-sky atmosphere (similar to the “opposition effect” in vegetation,e.g. Hapke et al. 1996). Figure 1b illustrates the processof entrapment by clouds when 3D transport is permitted:
radiation passing down through a clear-sky (or less optically thick) part of the atmosphere may be reflected backupward at a slantwise angle and encounter the base of acloud due to horizontal transport within either the clearsky or cloudy region. The depiction of “upward trapping”
by Várnai and Davies (1999) was similar except that thetwo cloud layers were part of the same cloud, and the reflected ray was intercepted by the edge rather than the baseof the upper layer. Since the area presented by the base of
a cloud is usually much larger than its edge, we would expect the impact of trapping by the base to be greater, on
average. Note that entrapment can also occur over reflective surfaces where the upward reflection is by the surface
rather than a cloud.
The interception of radiation by the upper cloud layer
reduces the reflectance of the scene, but the magnitudeof this effect depends on how far the radiation migrates
horizontally in the gap between the two cloud layers relative to the size of the clouds in the upper layer. Figure 1c
depicts the extreme case in which radiation is completely
horizontally homogenized in clear-sky layers. This “maximum entrapment” is actually the behavior of the original shortwave implementation of SPARTACUS describedby Hogan et al. (2016), as well as other solvers such as
the three-region solver in the original Edwards and Slingo(1996) radiation scheme that was adapted by Shonk and
Hogan (2008) to become the “Tripleclouds” solver. Shonkand Hogan (2008) described this radiative homogenization
as “anomalous horizontal transport,” which is not reallyaccurate as at least some of this transport occurs in reality.
Nonetheless, their method to remove it and thus to movefrom maximum entrapment (Fig. 1c) to zero entrapment(Fig. 1a) provides the starting point for representing more
realistic explicit entrapment (Fig. 1b) in SPARTACUS.
To demonstrate the importance of the difference between zero and maximum entrapment, we can make someidealizations and derive an analytic expression for scene

JOURNALOF

THE(a) Zero entrapment

3

ATMOSPHERIC SCIENCES(b) Explicit entrapment(c) Maximum entrapment

Layer 1
2
3
4
5
6
7

x 3.5

111111111
000000000
000000000
111111111

111111111
000000000
000000000
111111111

111111111
000000000
000000000
111111111
000000000
111111111

3.5

11111111
00000000
00000000
11111111
00000000
11111111

11111111
00000000
00000000
11111111
111111111
000000000
000000000
111111111
000000000
111111111

8

F IG . 1. Schematic illustrating SPARTACUS’s three possible treatments of “entrapment” of solar radiation underneath clouds in the case of two
randomly overlapped cloud layers each of cloud fraction 1/2. In each panel the degree to which the downwelling irradiance has been attenuatedby the clouds above is indicated by the darkness of the shading. The black-headed arrows depict representative light paths discussed in the text.
The double-headed arrow in panel b indicates the horizontal distance traveled by a single light ray reflected below half-level 3.5; the corresponding
mean horizontal distances are computed in section 4.

∗

2

2

R = R + T R/(1 − R ) = 2R/(1 + R).

(1)

Thus the scene reflectance in the zero-entrapment limit is
the weighted sum of the reflectance of the four columns:
Rzero = R/2 + R∗ /4 = R(1 + R/2)/(1 + R).

(2)

This is depicted by the solid line in Fig. 2. In the limit ofa perfectly reflective cloud (R = 1), the scene reflectancebecomes Rzero = 3/4.
Now consider the other extreme: maximum entrapment
as depicted in Fig. 1c. Since radiation passing througha clear-sky layer retains no memory of its horizontal location with respect to the clouds it has passed through inother layers, we consider the entire domain as a single column. Thus the reflectance of a layer with a cloud fractionof 1/2 is R/2, while its transmittance is 1− R/2. Applying
the Adding Method to obtain the reflectance of the scene
simply involves replacing R by R/2 in (1), yielding
Rmax = 2R/(2 + R).

(3)

This is depicted by the dashed line in Fig. 2, and in thelimit of a perfectly reflecting cloud the scene reflectancebecomes Rmax = 2/3. This is significantly less than the
ICA value of 3/4, suggesting that entrapment is an important process to treat when modeling the 3D shortwave
radiative effects of clouds.

0.8

Reflectance of scene

reflectance in both cases. As in Fig. 1, we consider thetwo cloud layers each to have a cloud fraction of 1/2, a
total cloud cover of 3/4 and to be in a vacuum over a surface with an albedo of zero. Both cloud layers have thesame reflectance R, and scatter conservatively so that theirtransmittance is T = 1 − R. In the zero-entrapment case
we apply the ICA, splitting the scene into four columns ofequal width. One is clear sky with a reflectance of zero,two contain a single cloud layer with reflectance R, andthe final column consists of two cloud layers, which the
Adding Method (Lacis and Hansen 1974) predicts to havea reflectance of

0.6

0.4

0.2
Zero entrapment (ICA)
Maximum entrapment
0

0

0.2

0.4

0.6

0.8

1

Reflectance of cloud layers
F IG . 2. Reflectance of the idealized scene discussed in section
2, composed of two randomly-overlapped non-absorbing cloud layers
each with a cloud fraction of 1/2, over a black surface in vacuum. The
solid line depicts Eq. 2 and the dashed line depicts Eq. 3.

3. Representing entrapment scenarios in SPARTACUS
Here we explain how SPARTACUS may be modified
to represent zero and maximum entrapment, illustrated inFig. 1, as well as the first step in representing explicit
entrapment. The symbols used in more than one equation in sections 3–5 are defined in Appendix B. SPARTACUS uses the Tripleclouds approach of splitting eachcloudy layer into three regions, one clear (denoted a) andtwo cloudy (denoted b and c) with different optical depths.
The radiation problem can then be written in terms of vectors and matrices; for example, u = (ua ub uc )T is a column vector containing the upwelling diffuse irradiances ata particular height in each of the three regions. Note thatthe irradiance component u j is defined to be power in region j per unit area of the entire grid box, not per unit areaof region j.

4

JOURNAL

OF

THE

ATMOSPHERIC SCIENCES

In order to represent the full range of entrapment scenarios, two aspects of the SPARTACUS implementation of
Hogan et al. (2016) require modification, both in the upward pass of the Adding Method (their Eqs. 24–30). The
first change is to describe the direct solar beam in termsof Di−1/2 , the albedo to direct radiation of the entire atmosphere and surface below half level i − 1/2. This accompanies Ai−1/2 , which is the corresponding albedo to
downwelling diffuse radiation. (As shown in Fig. 1a, we
index full atmospheric layers by i, counting down from thehighest layer i = 1, and half level i − 1/2 refers to the interface between layers i − 1 and i.) This change mirrors
the application by Hogan et al. (2018) of SPARTACUS tovegetation. Both of these albedos are matrices of the form
Aaa
A =  Aab
Aac


Aba
Abb
Abc


Aca
Acb  ,
Acc

(4)

where A jk is the fraction of diffuse downwelling radiationin region j that is reflected up in region k. These definitions ensure that u at any given height is equal to thesum of reflection of the downward diffuse irradiance v andthe downward direct irradiance s (Eq. 40 of Hogan et al.
2018): u = Av + Ds.
The second change needed to represent entrapment concerns how these albedo matrices are translated from the
top of one layer to the base of the layer above, accounting for the overlap of the clouds and associated regionsin the two layers. We follow Hogan et al. (2016) and define Abelow i−1/2 as the albedo of the atmosphere just below
half-level i − 1/2, so using the region definitions of layer
i. Likewise, Aabove i−1/2 is the albedo just above this halflevel, using the region definitions of layer i − 1. Equation
30 of Hogan et al. (2016) relates the two according to the
maximum-entrapment assumption:
Aabove i−1/2 = Ui−1/2 Abelow i−1/2 Vi−1/2 ,(5)where U and V are the upward and downward overlap
matrices. They are defined such that uabove = Uubelow
and vbelow = Vvabove ; that is, they map irradiances passing through a half-level on to the regions of the layer theradiation is entering. To see how (5) leads to maximum
entrapment, consider what happens at half-level 3.5 in the
two-region example of Fig. 1, where the matrices would
be:


1/2 0
U3.5 =
;
1/2 0


A 0
Abelow 3.5 =
;
0 0


1 1
V3.5 =,(6)
0 0

where A is the albedo of the atmosphere below half level
3.5. Applying (5) yields

A/2 A/2
Aabove 3.5 =.
(7)
A/2 A/2
This confirms that radiation exiting the base of either theclear or cloudy regions in layer 3 has an equal probabilityof being reflected back up into either of these two regions.
Shonk and Hogan (2008) described how to eliminate
this horizontal transport in their Tripleclouds solver andthereby achieve zero entrapment; their solution may be
written as
aabove i−1/2 = VTi−1/2 abelow i−1/2 ,(8)

where a is a column vector containing the reflectances ofeach region with the assumption that light is always reflected up from the same region it enters. Since Tripleclouds neglects lateral radiation flows between regions, A
is diagonal and a simply contains its diagonal elements.
To apply the zero-entrapment assumption to a SPARTACUS simulation that includes lateral flows between regions, abelow i−1/2 is defined such that its jth element contains the sum of the jth column of Abelow i−1/2 . Physicallythis means that flows represented by the white-headed arrows in Fig. 3, which involve reflection up into a different
region of the upper layer, are forced to be reflected up into
the same region. After applying (8), Aabove i−1/2 is defined
to be a diagonal matrix with its diagonal elements given
by aabove i−1/2 and the rest of SPARTACUS is unchanged.
Applying zero entrapment to the matrices in (6) yields

A 0Aabove 3.5 =
,
(9)
0 A
which indicates that, as required, radiation exiting the baseof either region of layer 3 is reflected back up into the same
region.
To represent explicit entrapment we first assume that radiation that travels far enough horizontally to pass througha cloud edge (represented by the off-diagonal elements ofAbelow i−1/2 and shown in Fig. 3a) has lost memory of theregion it originated from in layer i − 1, so can be treatedby maximum entrapment. By contrast, the destination ofreflected radiation that does not pass through a cloud edge
(the diagonal elements of Abelow i−1/2 illustrated in Fig. 3b)
is computed taking account of the mean horizontal distance traveled. Mathematically this is achieved (in the 2region case) by
Aabove i−1/2= Ui−1/2
+ ∑

j=a,b

0
Aab
below i−1/2jjj
Abelow i−1/2 Li−1/2Ababelow i−1/20!

aj
Vi−1/200Vi−1/2bj

Vi−1/2!

, (10)

JOURNALOFTHE(a) Off−diagonal elementsa5

ATMOSPHERIC SCIENCES(b) Diagonal elementsbaLayer:
b
i−1

i
b

a

b

b

a

b

F IG . 3. Schematic illustrating the radiation paths considered when computing entrapment at half-level i − 1/2. (a) The four paths that contributeto Aabbelowi−1/2 , which is the off-diagonal element of Abelow i−1/2 that represents the fraction of radiation downwelling into region a (clear sky) oflayer i that is scattered back up into region b (cloud) of the same layer. (b) The four paths that contribute to Aaa
below i−1/2 , which is the diagonal
element of Abelow i−1/2 that represents the fraction of radiation downwelling into region a of layer i that is scattered back up in the same region.
Note that in practice we include radiation paths that pass through the base of layer i and are reflected back by the lower layers. The white-headed
arrows represent radiation flows involving reflection back up into a different region of layer i − 1.

and similarly for Dabove i−1/2 . The first term on the righthand side is the maximum-entrapment formula (Eq. 5) applied just to the off-diagonal elements. The second termjj
on the right-hand side contains Abelow i−1/2 , the jth diagonal of Abelow i−1/2 , as well as elements from the jth row of
Vi−1/2 .
jMatrix Li−1/2 is central to SPARTACUS’s representation of explicit entrapment: it expresses how much lateral transport occurs for reflected radiation within region jof layer i, but also accounting for radiation passing down
through the layers below. Its elements quantify the weightof each of the arrows in Fig. 3b. Since it repartitions radiation between regions without changing the total energy,
its columns sum to 1 (i.e. it is a left stochastic matrix).
If we wished for the diagonal elements of Abelow i−1/2 tobe treated with zero entrapment (thereby eliminating thewhite-headed arrows in Fig. 3b), then we could simply setjLi−1/2 to the identity matrix leading to the second termon the right-hand side behaving exactly as (8). Otherwise,jLi−1/2 is computed in a two-step process. Section 4 describes and validates the first step: computing the mean
horizontal distance traveled by reflected radiation. Section
5 describes and validates the second step: using the mean
horizontal distance traveled to compute the elements ofjmatrix Li−1/2 and hence how much radiation is “trapped”
by passing across the dashed line in Fig. 3b.
4. Explicit entrapment: horizontal distance traveledby reflected radiation
The white-headed arrow in Fig. 1b illustrates the horizontal distance traveled by a single light ray reflected below half-level 3.5, and includes the horizontal distance associated with both the downward and upward parts of thejourney. This section deals with the task of estimating
the mean horizontal distance traveled by reflected radiation below a particular half-level i − 1/2, considering all
possible light paths, including those that penetrate down
and up through multiple layers. However, we exclude lightpaths that pass laterally through cloud boundaries, as these
are treated by maximum entrapment embodied in the firstterm on the right-hand side of (10). Since the downward
part of the journey depends on whether the incoming radiation at half-level i − 1/2 is diffuse or direct, we computeseparate diffuse and direct mean horizontal distances, denoted xi−1/2 and yi−1/2 , respectively. But note that lightrays contributing to yi−1/2 can only strictly be considered
to be “direct” sunlight until their first scattering event.
a. Method
We consider a plane-parallel atmosphere (i.e. containing one region in each layer), and later adapt it to multiple
regions. The horizontal distance traveled by direct radiation passing through a single layer i is∆yi = ∆zi tan θ0 ,(11)

where ∆zi is the physical thickness of the layer and θ0
is the solar zenith angle. It was shown by Schäfer et al.
(2016) that the equivalent expression for isotropic diffuse
radiation is∆xi = ∆zi π /2.
(12)
We use these as building blocks in our estimate of distancetraveled during reflection from multiple layers. Consider
how the two-stream Adding Method is applied for diffuse
radiation in a single layer i: given the layer reflectanceRi and transmittance Ti , the scene albedo at the top of thelayer is given by
Ai−1/2 = Ri +Ti2 Ai+1/2 1 + Ri Ai+1/2 + (Ri Ai+1/2 )2 + · · · ,(13)where Ai+1/2 is the scene albedo at the bottom of the layer
and the terms in the square brackets represent multiple reflections between layer i and the atmosphere below. Thisis a geometric series that reduces to(14)
Ai−1/2 = Ri + Ti2 Ai+1/2 / 1 − Ri Ai+1/2 .
It can be adapted to compute xi−1/2 , the mean horizontal
distance traveled by reflected radiation below the top of the

6

JOURNAL

OF

THE

ATMOSPHERIC SCIENCES

layer, as the weighted average of the distances associatedwith each order of scattering in (13):Ai−1/2 xi−1/2 = Ri x̂i + Ti2 Ai+1/2
× 2x̂i + xi+1/2h√i+ Ri Ai+1/2 x̂i + 2 x̂i + xi+1/2h√i+ (Ri Ai+1/2 )2 x̂i + 3 x̂i + xi+1/2+···},


Ti2 Ai+1/2 x̂i + xi+1/2
xi−1/2 ≃ x̂i +1.5 .
Ai−1/2 1 − Ri Ai+1/2(15)where xi+1/2 is the mean distance traveled by reflected radiation below the base of layer i. The first term on theright-hand side contains x̂i , the mean horizontal distancetraveled by radiation that is reflected by layer i, rather thanpenetrating the layer and being reflected by the layers below. We assume that, on average, such radiation penetrates to the center of the layer before being reflected backout (hence traveling a distance ∆xi /2 on each of the downward and upward paths), and that the azimuthal scattering
angle at the point of reflection is random (so these distances should be added in quadrature) leading to√1/2= ∆xi / 2.
x̂i = 2(∆xi /2)2

(16)This neglects the additional distance associated with multiple scattering entirely within a layer, but the good performance reported in section 4b suggests that this is a small
effect in practice.
The first line in the curly brackets in (15) represents
radiation that passes down through the entire layer i andback up again, so the horizontal distance associated withtransiting the layer is twice that of radiation reflected by
the layer (the 2x̂i term), and is added to the distance associated with reflection by the layers below half level i + 1/2(the xi+1/2 term). The subsequent lines in the curly brackets in (15) include multiple reflections between layer i andthe layers below half level i + 1/2. Each reflection addsx̂i + xi+1/2 to the distance traveled, but due to the random
azimuthal scattering angle, they again should be added inquadrature leading to the square-root weighting term.
Equation (15) may be rearranged to obtain
xi−1/2 = x̂i +
Ti2 Ai+1/2x̂i + xi+1/2
Ai−1/2∞ p
j× ∑ j + 1 Ri Ai+1/2 .

√∑∞j=0 j + 1 a j ≃ (1− a)−1.5 , which has errors of less than
10% for a < 0.9. Thus we approximate (17) by

(17)

j=0

Can we reduce this infinite series to a closed-form expression as before? Equation 14 exploited the fact that∑∞j=0 a j = (1 − a)−1 , and by differentiating this expression we find that ∑∞j=0 ( j + 1)a j = (1 − a)−2 . The infinite
series in (17) lies, in some sense, between these two series, and we find empirically that to a good approximation(18)

This equation may be applied sequentially from the surface up through the atmosphere to obtain a profile ofxi−1/2 . The surface value is xn+1/2 = 0, since reflection
right at the surface is not associated with any horizontaltransport.
Next we seek an equivalent expression for yi−1/2 , thehorizontal distance traveled by reflected direct radiation.
The equivalent expression to (14) for the albedo to direct
radiation is
S− Ai+1/2 + Ei Di+1/2,
Di−1/2 = Si+ + Ti i
1 − Ri Ai+1/2(19)where three new layer properties have been introduced: Ei
is the fraction of direct radiation that penetrates the layer
without being scattered, Si− is the fraction that penetrates
the layer but is scattered on the way (so emerging into the
diffuse stream at the base of the layer), while Si+ is thereflectance of the layer to direct radiation.
The mean horizontal transport associated with reflection by the Si+ term we denote as ŷi , and is assumed to be
associated with direct radiation that penetrates to the center of the layer followed by a scattering event and a diffuse
path back to the top of the layer. Therefore we add ∆yi /2and ∆xi /2 in quadrature:
ŷi =

1/2
1
∆y2i + ∆x2i
.
2

(20)

The two terms on the numerator of the right-hand side of(19) represent the two ways that direct radiation can penetrate the layer: with and without being scattered on theway. When scattering occurs, the diffuse albedo Ai+1/2 isinvolved, which is associated with diffuse horizontal distance xi+1/2 , whereas when scattering does not occur, thedirect albedo Di+1/2 is involved, which is associated withdirect horizontal distance yi+1/2 . Subsequent reflections,
governed by the denominator of (19), all involve the addition of the diffuse horizontal distance xi+1/2 . Applying
the logic of (15) leads to terms involving the addition, in
quadrature, of yi+1/2 (associated with a direct reflection)and one or more of xi+1/2 (associated with one or moreinternal diffuse reflections). Unfortunately, these terms do
not reduce conveniently to a closed-form expression like(18). Therefore we make the approximation that in any
term involving the combination of yi+1/2 and xi+1/2 , yi+1/2

JOURNAL

OF

THE

can be replaced by xi+1/2 , leading toyi−1/2 ≃ ŷi +

where

Ti

Di−1/2 −× Si Ai+1/2 ξ + Ei Di+1/2 (ξ − 1) x̂i + xi+1/2(21)+Ei Di+1/2 ŷi + yi+1/2 ,ξ = 1 − Ri Ai+1/2

−1.5

.

(22)

In a plane-parallel atmosphere, (14), (18), (19) and (21)may be applied sequentially from the surface to top-ofatmosphere (TOA), to obtain profiles of xi−1/2 and yi−1/2 .
We stress that even though the calculation proceeds in asingle upward pass through the atmosphere, the computed
mean horizontal distances include both the downward andupward parts of the journey of reflected light rays. We
need the horizontal distance traveled in partially cloudy
profiles, so seek vectors x and y whose jth elements contain the horizontal distances associated with region j. As
illustrated in Fig. 3b, each region is considered independently, so we may still use these four equations to step
the elements of x and y from the base of the layer to thetop. The other inputs to these equations, Ri , Ti , Ei , Si± ,Ai+1/2 and Di+1/2 are taken as the diagonal elements tothe corresponding matrices available in the SPARTACUS
computation. Physically, the diagonal elements are usedbecause we are interested in horizontal transport that remains within a region in layer i; radiation that passes laterally between regions in layer i was dealt with by (10).
The final aspect to deal with partially cloudy profiles is totranslate x and y from the regions below half-level i − 1/2
to the regions above. We use the relevant overlap matrix
similarly to the operation in (8):xabove i−1/2 = VTi−1/2 xbelow i−1/2 ,7

ATMOSPHERIC SCIENCES(23)

and likewise for y.
b. Evaluation
Here we evaluate the estimates of mean horizontal distance traveled by reflected radiation as a function of height(the values of x and y above), using Monte Carlo calculations by the model of Villefranque et al. (2019), which implements ray-tracing techniques from computer graphics
and permits the paths of individual photons to be tracked.
The results are shown in Fig. 4 for four cloud scenes
and three solar zenith angles in simulations using periodic
boundary conditions in the horizontal. All are at a single
wavelength in vacuum with idealized cloud optical properties over a Lambertian surface with an albedo of 0.2. The
first profile (panels a–d) consists of a plane-parallel cloudlayer containing isotropic scatterers with an optical depthof 1 and a single-scattering albedo of 0.999,999. Beneaththe cloud, all reflection is from the surface so the mean

horizontal distance traveled increases linearly with height
above the surface, z. The “direct” mean horizontal distance, y, increases with θ0 due to (11), but is not zero foroverhead sun since it includes the return journey from thesurface to height z as diffuse radiation. Within the cloud afraction of downwelling radiation is reflected by the cloud,rather than the surface, and so the mean horizontal distance is reduced. We see that the SPARTACUS estimatesusing the method described above are accurate to around
10% for x and 3% for y. The second profile (panels e–h) is
the same but with an optical depth of 5. The SPARTACUS
errors are somewhat larger at around 18%.
The last two profiles contain more realistic clouds. Bothassume an asymmetry factor of 0.86, appropriate for liquid clouds in the mid-visible. SPARTACUS then performs
the usual delta-Eddington scaling, treating some of theforward scattered light as if it had not been scattered atall. To achieve a fair comparison in terms of the definition of “direct” and “diffuse” radiation, the Monte Carlomodel takes the delta-Eddington-scaled extinction coefficient, and assumes a Henyey-Greenstein scattering phasefunction using the delta-Eddington-scaled asymmetry factor value of 0.462. Panels i–l of Fig. 4 show the results for a 6.4 × 6.4 km large-eddy simulation of cumulus clouds from Brown et al. (2002), which was also usedby Hogan et al. (2016) and is based on an observed casefrom the Atmospheric Radiation Measurement (ARM)
program. Panels m–p show the results for a 100 × 100 km
scene from a 250-m simulation by the Canadian Global
Environmental Multiscale (GEM) model of a multilayer
liquid cloud (Pacific scene 16). The GEM scenes are described in detail in section 6. In both the ARM and the
GEM cases, the typical SPARTACUS errors are 25% for x
and 6% for y. Given the simplifications involved, SPARTACUS performs very well in estimating horizontal distance traveled, and should be adequate to feed into the final step for computing entrapment.
5. Explicit entrapment:
trapped?

how much radiation isa. Method
Here we use the estimated mean horizontal distancesjfrom the previous section to compute the matrix Li−1/2 inj(10). The meaning of Li−1/2 in the two-region case canjbe explained by considering how it acts on wi−1/2 (x) =Tw ja w jb , where w jk is the fraction of radiative energy
in region j of layer i that lies beneath region k in layeri − 1. It is denoted as a function of the mean horizontal
distance traveled, x, since the radiation entered the layer.
jLi−1/2is then defined such that
j

j

j

j

wi−1/2 (xi−1/2 ) = Li−1/2 wi−1/2 (0),

(24)

8

JOURNAL

Slab optical depth 1
Height (km)4

Cloud fraction

4

(a)

4

SZA 44°

4

(c)

3

3

2

2

2

2

1

1

1

1

4

Slab optical depth 5
Height (km)SZA 0°
(b)

3

0
0

0.5

1

(e)

3

2

4

6

8
4

(f)

2

1

1

Direct Monte Carlo
Diffuse Monte Carlo
Direct SPARTACUS
Diffuse SPARTACUS

0
0

0.5

1
4

(i)

2

4

6

8

2

4

6

(g)

3

3

2

2

1

1

8

2

4

6

8

(k)

3

3

3

2

2

2

2

1

1

1

1

0
0

8

0.5

1
8

(m)

0
0

2

4

6

8

(n)

2

4

6

8

(o)

6

6

6

4

4

4

4

2

2

2

2

0
0

0.5

Cloud fraction

1

0
0

4

8

12

Horizontal distance (km)

8

2

4

6

8

2

4

6

8

(l)

0
8

6

0

6

0
0

8

4

(h)

0
4

3

0

2

0
0

4

(j)

0
4

0
0

SZA 70°
(d)

0
0

3

2

4

0
0

4

0

ARM Cumulus
Height (km)ATMOSPHERIC SCIENCES

THE

3

0

GEM Pacific
Height (km)

OF

(p)

0
0

4

8

12

0

Horizontal distance (km)

4

8

12

Horizontal distance (km)F IG . 4. The mean horizontal distances traveled by reflected light, for diffuse and direct downwelling radiation (x and y in the text, respectively),as computed from Monte Carlo simulations (solid lines) and estimated by SPARTACUS (dashed lines) using the method described in section 4.
Each row of panels corresponds to a different atmospheric profile described in the text, with the cloud fraction in the left column and subsequentcolumns showing results for solar zenith angles (SZAs) of 0◦ , 44◦ and 70◦ . The final row corresponds to the cloud shown in Fig. 5c. The error
bars show the error on the mean from the Monte Carlo calculations; note that the Monte Carlo model can only compute mean distances given a
sufficient number of photons, so does not show x above cloud top or y deep in an optically thick cloud.

where w j (0) represents the radiation partitioning at the instant the radiation enters the layer from above (when thejjhorizontal distance traveled is zero) and wi−1/2(xi−1/2) is
the radiation partitioning at the instant it leaves the layer
after being reflected (when the mean distance traveled isjjxi−1/2 ). Hence the diagonal elements of Li−1/2 are represented by the black-headed arrows in Fig. 3b and theoff-diagonals by the white-headed arrows. Note that (10)jjjincludes the albedo Abelow i−1/2 so Li−1/2 describes only

the redistribution of energy, and therefore the sum of thejelements in wi−1/2 (x) is constant at 1 for any x.
jHow can we compute Li−1/2 , accounting for the possibility of radiation passing across the dashed line in Fig.
3b more than once? This can occur via multiple scattering
events in the layers below i − 1/2, as well as via radiationpassing beneath small clouds and emerging on the otherside. The problem is analogous to the original SPARTACUS problem of working out the net exchange of radiation
through cloud sides, and so we can use the same method.

JOURNAL

OF

THE

ATMOSPHERIC SCIENCES

Essentially we wish to solve a system of coupled differential equations of the form:jdwi−1/2dx= Γ wi−1/2 ,j(25)

where the matrix Γ contains the rates of radiation exchange between the “subregions” of region j (with subregions defined by the regions of layer i − 1 that they arebeneath, illustrated by the dashed line in Fig. 3b), and− f ab + f ba
Γ=
,
(26)
+ f ab − f ba
where f kl is the rate at which radiation (direct or diffuse)is transported from subregion k to l, per unit increase inhorizontal distance traveled x. The solution to (25) is (24)jbut with Li−1/2 written as a matrix exponential:jj
Li−1/2 = exp Γ xi−1/2 .
(27)

The repeated elements in (26) mean that the matrix exponential may be computed efficiently for both the twoand three-region cases as described in the appendix of
Hogan et al. (2018). The same method is used to compute the contribution to the direct albedo matrix, but usingjthe direct horizontal distance traveled yi−1/2 .
b. Representing fractal behavior
To test the validity of this approach, we use the contrasting binary cloud scenes shown in Figs. 5a–5d, whichhave been generated by applying an optical-depth threshold to four of the GEM simulations described in section 6.
A scalar field is defined containing a value of 1 in the clear(black) areas and 0 in the cloudy (white) areas, which canbe thought of as solar radiation that has passed throughthe gaps between the clouds. Gaussian smoothing is thenapplied to the field with varying smoothing scales x, representing horizontal radiation transport beneath the cloud.
Previous studies of the interaction of radiation and clouds
have found a Gaussian to be reasonably good at describing
the horizontal distribution of diffuse radiation originatingfrom a point source (e.g. Hogan and Battaglia 2008; Wissmeier et al. 2013). The fraction of the total scalar field thatis then in the cloudy parts of the domain is the “trapped
fraction”, and is shown by the black lines in Figs. 5e–5h.
The dotted lines show the cloud cover, which corresponds
to the trapped fraction one would expect if the radiative
energy were completely homogenized horizontally (maximum entrapment).
To apply the matrix-exponential method describedabove to estimate the trapped fraction we need to define
the lateral exchange rates in (26). From the geometric arguments of Hogan and Shonk (2013), we would expect
f kl = Lkl /π ck ,

(28)

9

where Lkl is the length of the interface between regions
k and l (i.e. the perimeter length of the clouds) per unit
area of the domain, and ck is the fractional area of regionk, both of which can be obtained by analyzing the binary
cloud fields. We then apply (26) and (27) to obtain L, thebottom-left element of which is the fraction of radiative
energy beneath the cloud shown by the solid gray lines in
Figs. 5e–5h. It can be seen that this is a good prediction of
the trapped fraction for low values of x, but for larger x it
overestimates entrapment significantly.
This can be explained by the fractal nature of clouds.
The Hogan and Shonk (2013) definition of lateral exchange rates imposes a length scale on the cloud field
based on the perimeter length. For example, we could
define the “effective cloud scale” S (see Eq. 20 of
Hogan et al. 2018) such that normalized perimeter lengthisLab = 4ca (1 − ca )/S,

(29)

where S can be thought of as the size that equally-sizedsquares would need to have if they were placed randomly
on a grid and their fractional cover and total perimeter length were equal to the values for the actual cloudfield. The S values for the scenes in Fig. 5 are shown
above panels e–h. Substitution of (29) into (28) gives
f ab = 4(1 − ca )/π S, indicating that the rate of exchange
between regions is inversely proportional to S. Thus, if all
clouds indeed had a diameter of around S then we would
expect the trapped fraction to quickly approach the asymptotic value of the cloud cover for x ≫ S.
In reality the clouds span a wide range of scales, and thepresence of very large clouds reduces the trapped fraction
for larger x. Another way of looking at this is to recognize that since clouds are fractal, the effective perimeter
length, L̂ab , ought to be a power-law function of the lengthscale at which the cloud field is being probed. This is usually written as L̂ab ∝ x1−D , where D is the fractal dimension. Many studies have estimated the fractal dimensionof clouds, with D = 1.5 being a reasonable representative value (e.g. Cahalan and Joseph 1989; Gotoh and Fujii
1998; Wood and Field 2011), implying L̂ab ∝ x−1/2 .
The dashed gray lines in Fig. 5e–5h show the result of
making effective perimeter length depend on x in this way,with the formula leading to the best fit given by p
L̂ab = Lab × min 1, 0.4S/x .

(30)

This formula caps the effective perimeter length L̂ab to be
no larger than the measured value Lab for x < 0.4S. The
much better agreement with the curves computed from actual cloud fields gives us confidence that this formula is appropriate to use in computing entrapment within SPARTACUS. Note that we have not found it necessary to use (30)

10

JOURNAL(a) Atlantic case 10

Distance north (km)

100

THE

ATMOSPHERIC SCIENCES(b) Atlantic case 7100(c) Pacific case 16

100

75

75

75

50

50

50

50

25

25

25

25

0
0

25

50

75

100

0
0

Distance east (km)(e) S = 754 m

0.5

50

75

100

0
0

Distance east (km)(f) S = 1220 m

0.5

Gaussian
Simple model
Fractal model
Cloud cover limit

0.4

25

25

50

75

100

0

Distance east (km)(g) S = 2686 m

0.5

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.2

0.1

0.1

0.1

0.1

0

0
0

4

8

12

0
0

4

8

12

Horizontal distance traveled (km)

50

75

100

(h) S = 3842 m

0.5

0.4

Horizontal distance traveled (km)25

Distance east (km)

0.4

0.3

(d) Pacific case 25

100

75

0

Trapped fraction

OF

0
0

4

8

12

Horizontal distance traveled (km)

0

4

8

12

Horizontal distance traveled (km)F IG . 5. Evaluation of the part of the algorithm that computes the fraction of radiation “trapped” due to horizontal transport. (a–d) Binary cloudmasks obtained by applying a visible-optical-depth threshold to four GEM scenes. (e–h) The corresponding fraction of downwelling clear-sky
radiation that is trapped beneath a cloud as a function of the mean horizontal distance traveled, x; the black line shows the result of a Gaussiansmoothing of the cloud masks (treated as truth), while the gray lines show two candidate methods for SPARTACUS that take as input the cloudscale S (indicated above each of the lower panels) and cloud cover.

(a) Most overhangs: ζ =1
b

a
b

a

Layer:i−1i(b) Fewest overhangs: ζ =0
b

a

i−1

b

a

i

F IG . 6. Illustration of two model cloud configurations with the samecloud fraction and overlap parameter: (a) the configuration that maximizes the number of “overhangs” where explicit entrapment must be
calculated (shown by the dashed vertical lines with curved arrows acrossthem), and (b) the configuration that minimizes the number of overhangs. In SPARTACUS this is controlled by the ζ factor shown above
each panel.

also for computing lateral exchange across cloud boundaries, because the relevant length scale would be the horizontal distance traveled by radiation as it passes throughan individual model layer, which is much less than x.
c. Treatment of overhanging clouds
Applying these findings in SPARTACUS presents one
further issue to resolve, since as shown in Fig. 3b we are

not dealing with radiative exchange between regions, but
exchange between the “subregions” of region j in layeri, defined according to the regions above them in layeri − 1. Unfortunately the perimeter length of the interface
between these subregions is not completely defined by thevariables available to SPARTACUS. Consider the case of
two layers, each with a cloud fraction of 0.5, an overlap
parameter of αi−1/2 = 0.5 (i.e. halfway between maximum and random overlap) such that the combined cloud
cover of the two layers is 0.625, and a particular value forthe effective normalized cloud perimeter length in the top
layer L̂abi−1 . Figure 6a depicts the way in which the clouds
could be configured that maximizes the perimeter length
between the subregions of region a in layer i, which in this2D diagram is illustrated by the number of overhanging
clouds shown by the dashed lines.
Figure 6b depicts an alternative configuration at theother extreme: in the left half of the gridbox the clouds
are overlapped maximally and in the right half they areoverlapped randomly. Since the overhangs are only associated with the randomly overlapped half of the gridbox,
there are fewer overhangs. This treatment of cloud overlap
was explored by Shonk et al. (2010), who showed that thearea of the gridbox associated with maximum overlap ofjjregion j could be written as αi−1/2 × min(ci−1 , ci ). We areconcerned with removing the fraction of region j in layerji that is maximally overlapped, so divide through by ci toJOURNAL

OF THEATMOSPHERIC SCIENCES

obtain the fraction of region j that is randomly overlapped:
Ci−1/2 = 1 − αi−1/2 × min(ci−1 , ci )/ci
j

j

j

j

(31)

Thus the most overhangs (Fig. 6a) is obtained by usingL̂kl
i−1 for the effective perimeter length between regions k
and l, while the fewest overhangs (Fig. 6b) is obtained byjusing L̂kl
i−1 × Ci−1/2 . There is no theoretical or observational evidence to support which is the most likely, so we
introduce a user-defined “overhang factor” ζ that varies
the effective perimeter length linearly between most overhangs (ζ = 1) and fewest overhangs (ζ = 0). This factor isleft as a parameter to be determined in section 6 according
to which leads to SPARTACUS predicting 3D radiative effects most accurately. It turns out that this property is ofsecond-order importance compared to whether entrapment
is treated at all.
6. Results
In this section we evaluate the shortwave 3D radiative effect predicted by the new SPARTACUS implementation in the ecRad radiation scheme (Hogan and Bozzo2018), and investigate the impact of various different treatments of entrapment. We have used 65 scenes generated
from simulations by Environment and Climate Change
Canada’s Global Environmental Multiscale model (GEM,
Girard et al. 2014), using the configuration described by
Leroyer et al. (2014) with the Milbrandt et al. (2005)
double-moment bulk microphysics cloud scheme.
Each scene measures 100×100 km, has a horizontal resolution of 250 m and employs 56 vertical layers. The simulations were originally performed to generate synthetic
satellite data from two swaths: an Atlantic swath on 7 December 2014 from Greenland to the Dominican Republic,
from which 39 scenes were extracted, and a Pacific swathon 24 June 2015 from Hawaii to Tonga, from which a further 26 scenes have been extracted. Thus the scenes spana wide range of cloud conditions.
Both ICA and 3D Monte Carlo shortwave radiativetransfer calculations have been performed on these scenes
using the model of Barker et al. (2003), which tracks photons through sequences of scattering events until they are
either absorbed by a particle, molecule, or the surface, orexit the domain’s top. Calculations were performed forsolar zenith angles at 5◦ intervals between 0◦ and 85◦ , butwith random solar azimuth angle, and assuming a periodic
domain. The Rapid Radiative Transfer Model for GCMs
(RRTM-G) of Iacono et al. (2008) was used to represent
gas absorption, the Yi et al. (2013) scheme for ice optical properties and Mie theory for liquid droplets. Scattering by air molecules and cloud particles were handled by
the Rayleigh and Henyey-Greenstein phase functions, respectively. To simplify the comparison with 1D radiation
schemes, all calculations assumed a Lambertian surface
with an albedo of 0.05.

11

We first compare the TOA cloud radiative effect between Monte Carlo calculations run in an ICA mode
(Monte Carlo ICA) and Tripleclouds (the SPARTACUS
control) for the same scenes, i.e. in the absence of 3D radiative transfer. In addition to cloud fraction and gridboxmean liquid and ice mixing ratio, Tripleclouds takes as
input the overlap parameter at each half-level and thefractional standard deviation of in-cloud water content,
FSD, in each layer. It was found that the original implementation of Tripleclouds was not capable of accurately
representing the effect of horizontal heterogeneities for
FSD > 1.5, which occurs in many of these scenes. Appendix A describes an improvement to Tripleclouds and
SPARTACUS that has overcome this problem. Figures
7a–7c reveal that the resulting root-mean-squared error in
CRE predicted by Tripleclouds is around 10% and its bias
is only 1–2%.
The differences between 3D and ICA Monte Carlo calculations of CRE for the 65 scenes are summarized by theblack box-and-whisker plots in Figs. 8a and 8b, and themean by the thick black line. We see that for θ0 in therange 0–75◦ , 3D effects increase both TOA and surface
CRE (i.e. make them less negative), and therefore act inthe sense of warming the Earth system by making clouds
less reflective. In individual cases, 3D effects can act tocool the Earth system by up to 7 W m−2 at θ0 = 75◦ , viainterception of direct solar radiation by cloud sides, but
plenty of other scenes have a warming effect even for large
θ0 .
The various SPARTACUS simulations, which all use
Tripleclouds as their 1D control, enable us to elucidate
the role of entrapment in explaining this behavior. SPARTACUS has been run taking as input the observed cloud
perimeter length in each layer, Lab , and the length of thecontour separating the optically thinner and thicker partsof the in-cloud extinction distribution, Lbc . In practice theformer is characterized by the cloud effective scale definedby (29), and the latter by a cloud heterogeneity scale, Shet ,given analogously by
Lbc = 4cc (1 − cc )/Shet .

(32)

The blue lines in Figs. 8a and 8b show the zeroentrapment scenario, which was illustrated in Fig. 1a. The
3D effect is much weaker overall; at TOA, the greatest
mean cooling is 1.9 W m−2 at θ0 = 75◦ , and the greatest mean warming is 1.7 W m−2 at θ0 = 0◦ . This general pattern can be explained by the mechanisms of shortwave side illumination and downward escape discussed by
Hogan and Shonk (2013), but downward escape is clearly
too weak a mechanism to explain the strong 3D effectfound in the Monte Carlo simulations for overhead sun.
At the other extreme, the green lines show the maximum
entrapment scenario, which was illustrated in Fig. 1c andinvolves complete horizontal homogenization of radiation

OF

Solar zenith angle: 20°
(a)

-200

-400

-600
-2
-2

RMSE: 29.0 W m-2

-800

-600

-400

-200

(b)
-200

-400
-2

-600

RMSE: 20.5 W m-2

0

-600

-2

)
(d)

-150

-2
-2

Bias: +1.2 W m
RMSE: 5.2 W m -2

-150

-100

-50

0
-2

Monte Carlo ICA CRE (W m )

(e)

20
(f)

30

40

10

30

20

20

10
0

10
MC mean: +6.7 W m -2

0
-10
-10

-100

-200
-200

0

40

-2

)

-200

-50

Monte Carlo ICA CRE (W m )

60

-2

-400

(c)

-2

Monte Carlo ICA CRE (W m )50-2Solar zenith angle: 80°

0

)

-800

Solar zenith angle: 50°

0

-2

0

ATMOSPHERIC SCIENCES

THE

Tripleclouds 1D CRE (W m -2)

JOURNAL

Tripleclouds 1D CRE (W m -2)

Tripleclouds 1D CRE (W m -2)

120

MC mean: +4.3 W m -2

-2

-2

Bias: +0.3 W m
RMSE: 1.5 W m -2

-2

Bias: +0.3 W m
RMSE: 6.8 W m -2

0

10

20

30

40

50

-10
-10

60
-2

-2

Bias: +0.3 W m
RMSE: 4.3 W m -2

0

10

20

30
-2

)

-10
-10

40

0

10

20
-2

)

)

F IG . 7. (a–c) Comparison of top-of-atmosphere shortwave cloud radiative effect (CRE) between reference Independent Column Approximation(ICA) calculations using the Monte Carlo (MC) model on the horizontal axis, and Tripleclouds (the SPARTACUS control) on the vertical axis, foreach of the 65 GEM scenes at three different solar zenith angles indicated above the panels. (d–f) Corresponding comparison of the 3D radiative
effect, i.e. the difference in CRE between 3D and 1D calculations. The SPARTACUS calculations use explicit entrapment with an overhang factorof ζ = 0, also shown by the red lines in Fig. 8. Each panel also states the mean of the MC calculations, and the bias and root-mean-squared error(RMSE) of the SPARTACUS/Tripleclouds calculations with respect to MC.
(a) Top-of-atmosphere CRE40(b) Surface CRE(c) Atmospheric absorption
Monte Carlo
SPARTACUS Zero
SPARTACUS Explicit (ζ=0)
SPARTACUS Explicit (ζ=1)
SPARTACUS Maximum

35

3D change (W m -2)

30
25
20
15
10
5
0
-5
-10
0

15

30

45

60

Solar zenith angle (°)

75

90

0

15

30

45

60

Solar zenith angle (°)

75

90

0

15

30

45

60

75

90

Solar zenith angle (°)

F IG . 8. Difference between simulations with and without 3D effects in (a) top-of-atmosphere cloud radiative effect (CRE), (b) surface CRE and(c) atmospheric absorption. The thick black and red lines show the average of the 65 GEM scenes for the Monte Carlo model and SPARTACUS in
its preferred configuration of explicit entrapment with ζ = 0. The box plots represent the corresponding median and interquartile range of the 65
scenes, with the whiskers representing the 5th and 95th percentiles. The other thick lines show averages of different SPARTACUS configurations.

JOURNAL

OF THEATMOSPHERIC SCIENCES

in each clear or cloudy region. This has a strong warming effect, reaching 19.6 W m−2 at TOA for overhead sun,but is over twice as strong as the reference Monte Carlocalculations.
The other two SPARTACUS simulations in Figs. 8a and8b are much closer to the reference calculations: the red
and pink lines show results using explicit entrapment described in sections 4 and 5, with the two treatments ofcloud overhangs illustrated in Fig. 6. It is clear that theleast-overhang scenario (ζ = 0) agrees with the Monte
Carlo calculations best at TOA; there is still a slight overestimate of the warming effect of 3D radiative transfer,but it is less than 8% for θ0 ≤ 55◦ . Therefore, ζ = 0 is ourpreferred SPARTACUS configuration for the remainder ofthe paper. Figures 7d–7f compare individual 3D radiative
effects between Monte Carlo and this SPARTACUS configuration, and while there is some scatter, the correlation
coefficients of 0.79–0.86 confirm that there is skill in predicting 3D effects for individual cases.
Figure 8c shows the change to total atmospheric absorption when 3D effects are included. The Monte Carlo calculations show an increase in absorption by around 1 W m−2
at most solar zenith angles, which is comparable to thefindings of Barker et al. (2016). Both the main 3D mechanisms contribute to this effect: side illumination at large
solar zenith angle enhances the interception and hence absorption by clouds, while Fig. 3b shows that entrapment
increases the path length of radiation in clear skies beneath cloud, enhancing water-vapor absorption. SPARTACUS with explicit entrapment leads to around 2 W m−2
greater atmospheric absorption than Tripleclouds, on average, which is twice the 3D effect in the Monte Carlosimulations. This is related to the presence of a handfulof outliers amongst the SPARTACUS simulations (shownby red dots in Fig. 8c); indeed, if we were to look at themedian rather than the mean of the 65 cases then it would
suggest instead that SPARTACUS tends to underestimate
the 3D effect on atmospheric absorption.
To investigate the factors that influence the nature of 3D
radiative transfer in individual cases, and the fidelity withwhich they are captured by SPARTACUS, we analyze theradiation fields for the four contrasting GEM scenes depicted in Fig. 9. Vertical profiles of the four main inputs to
SPARTACUS are shown in Figs. 10a–10d. Atlantic case
6 consists of cumulus clouds with some vertical development; the small effective cloud scale of S ≃ 1 km, andhence large cloud-side area, leads to significant shortwave
side illumination, with Fig. 10e showing a 7 W m−2 increase in the reflectance of the scene at θ0 = 75◦ (a cooling
effect). The dependence on θ0 is well captured by SPARTACUS, including the change in sign to a 7 W m−2 warming at θ0 = 0◦ , or a −9% change to CRE (see the captionof Fig. 10 for the total CRE values for overhead sun), but
the similarity between the various entrapment configura-

13

tions highlights that in this case the 3D effect for high sun
is mainly due to downward escape.
The three remaining scenes, by contrast, appear to be
dominated by entrapment. Atlantic case 14 contains deep
frontal cloud with considerable small-scale structure. The
zero-entrapment simulation in Fig. 10f shows the significant cooling effect of side illumination, but the explicitentrapment simulation shows that this is overwhelmed byentrapment and indeed the net warming by 3D effects is
up to 39 W m−2 (a −8% change to CRE). A key factor
is the large vertical extent of the cloud, which means thatradiation passing down through the gaps in the clouds cantravel a large distance horizontally before being reflectedback up to its original level, increasing the trapping. Atlantic case 32 contains much more homogeneous and overcast boundary-layer cloud. The zero-entrapment simulation has a 3D radiative effect of less than 1 W m−2 , confirming that cloud-side effects are weak. With entrapment
included, the 3D effect is a warming of up to 6–7 W m−2 ,with good agreement between SPARTACUS and Monte
Carlo. In absolute terms this effect is significant, but thisscene is the most reflective of the four and in relative terms
it is only a −1% change to CRE. Additional Tripleclouds
and SPARTACUS calculations in which the in-cloud heterogeneity is removed (i.e. setting FSD = 0) lead to the3D effect almost entirely disappearing, which suggest that
it is due to trapping associated with cloud heterogeneity,
similar to one of the mechanisms proposed by Várnai andDavies (1999). Finally, Pacific case 25 consists of remnants of deep convection including anvils with S ≃ 10 km.
Again the entrapment mechanism appears to dominate.
7. Conclusions
Cloud scenes have varied and complex structures, andconsequently it can be very challenging to understand themagnitude and even sign of the differences between radiation calculations with and without horizontal transport.
The simplest mechanism to understand shortwave 3D radiative transfer is side illumination, which enhances cloudreflectance. This has led many previous studies to focus
on cloud types with a relatively large cloud-side area suchas cumulus (Benner and Evans 2001; Pincus et al. 2005)
and aircraft contrails (Gounou and Hogan 2007). However, Barker et al. (2015) analyzed a much more varied andrepresentative set of cloud fields and found that shortwave
3D transport tends to reduce the reflectance of clouds overall and hence has a warming effect on the Earth system. In
this paper we propose the mechanism of “entrapment” toexplain this behavior. Entrapment is similar to one of themechanisms suggested by Várnai and Davies (1999) for
single-layer cloud scenes, but an important insight is thatit need not involve transport through cloud sides. It tends
to be strongest in deep, multilayer scenes, which are common in reality but have tended to be ignored in previous

14

JOURNAL

OF

THE

ATMOSPHERIC SCIENCES

F IG . 9. Three-dimensional visualizations of four contrasting GEM scenes, where the grid axes are marked in kilometers. The results of 3D radiation
calculations on these scenes are shown in Fig. 10.
Atlantic case 6
Fractional standard deviation

Atlantic case 14
Fractional standard deviation

Atlantic case 32
Fractional standard deviation

Pacific case 25
Fractional standard deviation

0

0

0

0

1

2

3

4

1

2

3

4

1

2

3

4

1

2

3

4

15

Height (km)

(a)

(b)

(c)

(d)

10
Frac. standard dev.
Cloud scale
Heterogeneity scale

5

0
0

10

20

30

0

Cloud scale (km)

10

20

30

0

Cloud scale (km)

10

20

30

0

Cloud scale (km)

10

20

30

Cloud scale (km)3D change to CRE (W m -2)

50
(e)

40

(f)

(g)

(h)

SPARTACUS Zero
SPART. Explicit (ζ=0)
SPART. Explicit (ζ=1)
SPART. Maximum
Monte Carlo

30
20
10
0
-10
0

30

60

Solar zenith angle (°)

90

0

30

60

Solar zenith angle (°)

90

0

30

60

Solar zenith angle (°)

90

0

30

60

90

Solar zenith angle (°)

F IG . 10. (a–d) The four main cloud geometry parameters used by SPARTACUS for the four GEM scenes shown in Fig. 9, where cloud fraction ingray varies between 0 and 1 across the horizontal axis, and the cloud and heterogeneity scales are defined by (29) and (32). (e–h) The corresponding
calculations of the difference in top-of-atmosphere shortwave cloud radiative effect (CRE) between 3D and 1D calculations, where the black lines
show Monte Carlo results and the colored lines show various configurations of SPARTACUS. For reference, the overhead-sun CRE for the four
scenes according to the 3D Monte Carlo model was −79, −476, −832 and −184 W m−2 , respectively.

case studies, presumably due to them being regarded as
too complex to interpret.

We have described modifications to the shortwave
SPARTACUS solver of Hogan et al. (2016) to incorporate

OF THE

15

ATMOSPHERIC SCIENCES

an explicit calculation of entrapment, making use of theeffective cloud scale variable already provided as input to
SPARTACUS. This involves a novel method to estimate
the mean horizontal distance traveled by reflected radiation, something that could be useful in other contexts, forexample in determining when the radiation scheme of acloud-resolving model ought to represent lateral exchangeof radiation between gridboxes. We have also found it
necessary to explicitly represent the fractal dimension ofcloud perimeters.
Evaluation against Monte Carlo calculations on 65 contrasting scenes from a cloud-resolving model reveals the
new SPARTACUS scheme to be capable of predicting
the “3D effect”, i.e. the difference between cloud radiative effect computed with and without horizontal radiativetransport, with a TOA bias of no more than 0.3 W m−2
for all solar zenith angles, and skill in predicting the dependence of the 3D effect on solar zenith angle in individual scenes. On average, 3D radiative effects tend to
make these scenes less reflective (similar to the findings ofBarker et al. 2015), implying that entrapment is a more important mechanism than side illumination. However, thisresult is highly dependent on the realism of the clouds simulated by the cloud-resolving model; if real clouds were
smaller, on average, than those used here then the sideillumination mechanism would be relatively more important.
The modified SPARTACUS is now an option in theecRad radiation scheme (Hogan and Bozzo 2018) used in
the ECMWF model. Hogan and Bozzo (2018) reported
the original SPARTACUS with maximum entrapment tobe 3.5 times slower than Tripleclouds, and we find that explicit entrapment increases this to around 4.5. While too
costly to use operationally, it is fast enough to use for research purposes. The next step will be to use this validated
tool to estimate the global impact of 3D radiative transfer,
not just in the shortwave but also in the longwave (Schäfer
et al. 2016).
Acknowledgments
We thank Zhipeng Qu for performing the GEM model
simulations. NV acknowledges support from the Agence
Nationale de la Recherche (grant ANR-16-CE01-0010),
and from the French Ministry of Higher Education, Research and Innovation through the doctoral school SDU2E
of Université de Toulouse.
Appendix A: Improving Tripleclouds for very heterogeneous scenes
SPARTACUS can be thought of as a 3D extension of the
1D Tripleclouds solver of Shonk and Hogan (2008), so
before evaluating SPARTACUS against 3D Monte Carlocalculations we need to be sure that Tripleclouds agrees

Fractional standard deviation

JOURNAL

4

4

3

3

2

2

1

1
(a)

0

0

(b)
20

0

40

0

20

Mean optical depth
-0.1

40

Mean optical depth

-0.05

0

0.05

0.1

Albedo bias

F IG . 11. Albedo bias of the Tripleclouds method for a gamma distribution of cloud optical depths with increasing fractional standard deviation indicated on the vertical axis: (a) equal-area cloudy regions withthe optical depth of the first cloudy region given by the 16th percentileof the gamma distribution; (b) new method in which the area of the two
cloudy regions can be different. Following Shonk and Hogan (2008),
this analysis approximates the relationship between albedo, α , and optical depth, δ , as α = 0.2 + 0.525δ /(δ + 3.5).

well with the corresponding ICA calculations. Tripleclouds takes as input a profile of “fractional standard deviation” (FSD), which is the standard deviation of in-cloudwater content or extinction coefficient divided by the incloud mean value. Tripleclouds divides the cloud in eachlayer into two regions (denoted b and c) of different extinction coefficient. Shonk and Hogan (2008) reported thatfor FSD up to 2, predicted irradiances agreed best withICA if the two cloudy regions had equal area, region b
used the 16th percentile of the full extinction distribution,and the extinction of region c was chosen so as to conserve the layer-mean extinction. The implementation of
Tripleclouds in ecRad (Hogan and Bozzo 2018) includes
the option to represent either a lognormal distribution of
optical depth, in which case the ratio of the 16th percentile to the in-cloud mean is given approximately by Eq.
44 of Hogan et al. (2016), or a gamma distribution (e.g.
Barker et al. 1996) for which this ratio is approximated by
rb = exp −FSD − FSD2 /2 − FSD3 /4 .

(33)

Conservation of mean extinction coefficient requires theratio of the extinction of region c to the in-cloud mean to
be rc = 2 − rb .
Comparison of Tripleclouds to ICA calculations on thescenes described in section 6 revealed the gamma distribution to perform best, but even then Tripleclouds tended tooverestimate scene reflectance for the more heterogeneous
scenes, some of which have FSD values up to 4. In Fig.
11a we have repeated the analysis of Shonk and Hogan(2008) but for a gamma rather than a lognormal distribution, and considered larger values of FSD. A substantial

16

JOURNAL

OF

THE

ATMOSPHERIC SCIENCES

albedo overestimate is apparent for FSD > 2. The problem arises because for large FSD, rb tends to 0 and rc to
2. Since the two cloudy regions have equal area, the actual fractional standard deviation of the Tripleclouds representation of the in-cloud extinction distribution tends to
1, which may be much less than the FSD of the gamma
distribution being approximated.
The solution we propose to overcome this problem is
twofold. First, a new rb is defined with a lower limit of
0.025:brnew= 0.975rb + 0.025.
(34)
Second, for large values of FSD we increase the fractional
area of region b and correspondingly reduce that of region
c: for FSD in the range 1.5–3.75, the fraction of the cloudy
area occupied by region b increases linearly from 0.5 to
0.9, while outside this range it is capped at 0.5 or 0.9. The
extinction of region c is still chosen to conserve the layermean value. Figure 11b shows that these changes virtuallyeliminate the albedo bias up to an FSD of 4. This solution has been implemented in both the Tripleclouds and
SPARTACUS solvers of ecRad, and is used in section 6.
Appendix B: List of symbols
The following list includes symbols used in more than
one equation in sections 3–5, and “PP” indicates a variable from section 4a where a plane-parallel atmosphere
has been assumed.
Aabove i−1/2 diffuse albedo of entire atmosphere and surface below interface i − 1/2, with matrix elements configured for regions in the layer
above the interface (region i − 1)
Ai−1/2 diffuse albedo of entire atmosphere and surface below interface i − 1/2 (PP)ck fraction of layer occupied by region k

Dabove i−1/2 as Aabove i−1/2 but for direct radiation
Di−1/2 direct albedo of entire atmosphere and surface below interface i − 1/2 (PP)

Ri diffuse reflectance of layer i (PP)s vector of downwelling direct irradiances ineach region at a particular heightSi− reflectance of layer i to direct radiation (PP)Si+ fraction of direct radiation that penetrates
layer i and is scattered on the way (PP)S, Shet effective cloud scale; cloud heterogeneityscale
Ti transmittance of layer i to diffuse radiation(PP)u vector of upwelling irradiances in each region at a particular height
Ui−1/2 upward overlap matrix expressing how upwelling irradiances in each region just below
interface i − 1/2 are transported into the regions just above
v vector of downwelling diffuse irradiances ineach region at a particular height
Vi−1/2 downward overlap matrix expressing how
downwelling irradiances in each region justabove interface i − 1/2 are transported into
the regions just belowjwi−1/2 vector expressing the fraction of radiation inregion j of layer i that is beneath each regionof layer i − 1
xi−1/2 mean horizontal distance traveled by reflected diffuse radiation below interface i −
1/2 (PP)x̂i mean horizontal distance traveled by diffuse
radiation reflected by layer i (PP)yi−1/2 mean horizontal distance traveled by reflected direct radiation below interface i −
1/2 (PP)Ei fraction of direct radiation penetrating layer
i without being scattered (PP)ŷi mean horizontal distance traveled by direct
radiation reflected by layer i (PP)f kl rate at which radiation passes from subregion k to l, per unit increase in horizontal
distance traveledΓ matrix expressing the rates of radiation exchange between the subregions of a region

Lkl length of interface between regions k and lnormalized by the area of the domain∆xi horizontal distance traveled by diffuse radiation passing through layer iL̂kl effective normalized interface length, accounting for the fractal nature of clouds∆yi horizontal distance traveled by direct radiation passing through layer ijLi−1/2matrix expressing how much radiation entering region j of layer i from the top and subsequently reflected back up, is exchanged
between the various regions in layer i − 1∆zi thickness of layer iζ overhang factor
θ0 solar zenith angle

JOURNAL

OF THEATMOSPHERIC SCIENCES

References
Barker, H. W., and J. A. Davies, 1992: Solar radiative fluxes for broken
cloud fields above reflecting surfaces. J. Atmos. Sci., 49, 749–761.
Barker, H. W., B. A. Wielicki and L. Parker, 1996: A parameterization
for computing grid-averaged solar fluxes for inhomogeneous marine
boundary layer clouds: 3. Validation using satellite data. J. Atmos.
Sci., 53, 2304–2316.
Barker, H. W., R. K. Goldstein and D. E. Stevens, 2003: Monte Carlosimulation of solar reflectances for cloudy atmospheres. J. Atmos.
Sci., 60, 1881–1894.
Barker, H. W., J. N. S. Cole, J. Li, B. Yi and P. Yang, 2015: Estimationof errors in two-stream approximations of the solar radiative transfer
equation for cloudy-sky conditions. J. Atmos. Sci., 72, 4053–4074.
Barker, H. W., J. N. S. Cole, J. Li and K. von Salzen, 2016: A
parametrization of 3-D subgrid-scale clouds for conventional GCMs:
Assessment using A-Train satellite data and solar radiative transfer
characteristics, J. Adv. Model. Earth Syst., 8, 566–597.
Benner, T. C., and K. F. Evans, 2001: Three-dimensional solar radiativetransfer in small tropical cumulus fields derived from high-resolution
imagery. J. Geophys. Res., 106, 14,975–14,984.
Brown, A. R., R. T. Cederwall, A. Chlond, P. G. Duynkerke, J. Golaz,
M. Khairoutdinov, D. C. Lewellen, A. P. Lock, M. K. MacVean, C.
Moeng, R. A. Neggers, A. P. Siebesma and B. Stevens, 2002: Largeeddy simulation of the diurnal cycle of shallow cumulus convection
over land. Q. J. R. Meteorol. Soc., 128, 1075–1093.
Cahalan, R. F., and J. H. Joseph, 1989: Fractal statistics of cloud fields.
Mon. Wea. Rev., 117, 261–272.
Di Giuseppe, F., and A. M. Tompkins, 2015: Generalizing cloud overlap treatment to include the effect of wind shear. J. Atmos. Sci., 72,
2865–2876.
Edwards, J. M., and A. Slingo, 1996: Studies with a flexible new radiation code: 1. Choosing a configuration for a large scale model. Q. J.
R. Meteorol. Soc., 122 689–719
Girard, C., M. Desgagné, R. McTaggart-Cowan, J. Côté, M. Charron,
S. Gravel, V. Lee, A. Patoine, A. Qaddouri, M. Roch, L. Spacek,
M. Tanguay, P. A. Vaillancourt and A. Zadra, 2014: Staggered vertical discretization of the Canadian Environmental Multiscale (GEM)model using a coordinate of the log-hydrostatic-pressure type. Mon.
Wea. Rev., 142, 1183–1196.
Gotoh, K., and Y. Fujii, 1998: A fractal dimensional analysis on thecloud shape parameters of cumulus over land. J. Appl. Meteorol., 37,
1283–1292.
Gounou, A., and R. J. Hogan, 2007: A sensitivity study of the effect ofhorizontal photon transport on the radiative forcing of contrails. J.
Atmos. Sci., 64, 1706–1716.
Hapke, B., D. DiMucci, R. Nelson and W. Smythe, 1996: The causeof the hot spot in vegetation canopies in soils: shadow-hiding versus
coherent backscatter. Rem. Sens. Env., 58, 63–68.
Hill, P. G., C. J. Morcrette and I. A. Boutle, 2015: A regime-dependent
parametrization of subgrid-scale cloud water content variability. Q.
J. R. Meteorol. Soc., 141, 1975–1986.
Hogan, R. J., and A. Battaglia, 2008: Fast lidar and radar multiplescattering models: 2. Wide-angle scattering using the timedependent two-stream approximation. J. Atmos. Sci., 65, 3636–3651.
Hogan, R. J., and J. K. P. Shonk, 2013: Incorporating the effects of 3D
radiative transfer in the presence of clouds into two-stream multilayer radiation schemes. J. Atmos. Sci., 70, 708–724.
Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation
scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10,
1990–2008.
Hogan, R. J., S. A. K. Schäfer, C. Klinger, J.-C. Chiu and B.
Mayer, 2016: Representing 3D cloud-radiation effects in two-stream
schemes: 2. Matrix formulation and broadband evaluation. J. Geophys. Res. Atmos., 121, 8583–8599.

17

Hogan, R. J., T. Quaife and R. Braghiere, 2018: Fast matrix treatment of 3-D radiative transfer in vegetation canopies: SPARTACUSVegetation 1.1. Geosci. Model Devel., 11, 339–350.
Iacono, M. J., J. S. Delamere, E. J. Mlawer, M. W. Shephard, S. A.
Clough and W. D. Collins, 2008: Radiative forcing by long-lived
greenhouse gases: Calculations with the AER radiative transfer models. J. Geophys. Res., 113, D13103, doi:10.1029/2008JD009944.
Lacis, A. A., and J. E. Hansen, 1974: A parameterization for the absorption of solar radiation in the Earth’s atmosphere. J. Atmos. Sci.,
31, 118–133.
Leroyer, S., S. Bélair, S. Husain and J. Mailhot, 2014: Subkilometer
numerical weather prediction in an urban coastal area: A case study
over the Vancouver metropolitan area. J. Appl. Meteor. Climatol., 53,
1433–1453.
Milbrandt, J. A., and M. K. Yau, 2005: A multi-moment bulk microphysics parameterization: 3. A proposed three-moment closure andscheme description. J. Atmos. Sci., 62, 3065–3081.
Pincus, R., H. W. Barker and J.-J. Morcrette, 2003: A fast,
flexible, approximate technique for computing radiative transfer in inhomogeneous cloud fields. J. Geophys. Res., 103,
doi:10.1029/2002JD003322.
Pincus, R., C. Hannay and K. F. Evans, 2005: The accuracy of determining three-dimensional radiative transfer effects in cumulus clouds
using ground-based profiling instruments. J. Atmos. Sci., 62, 2284–
2293.
Schäfer, S. A. K., R. J. Hogan, C. Klinger and B. Mayer, 2016: Representing 3D effects in two-stream radiation schemes: 1. Longwave
considerations and effective cloud edge length. J. Geophys. Res. Atmos., 121, 8567–8582.
Shonk, J. K. P., and R. J. Hogan, 2008: Tripleclouds: An efficientmethod for representing horizontal cloud inhomogeneity in 1D radiation schemes by using three regions at each height. J. Climate,21, 2352–2370.
Shonk, J. K. P., R. J. Hogan, J. M. Edwards and G. G. Mace, 2010:
Effect of improving representation of horizontal and vertical cloud
structure on the Earth’s global radiation budget: 1. Review andparametrization. Q. J. R. Meteorol. Soc., 136, 1191–1204.
Várnai, T., and R. Davies, 1999: Effects of cloud heterogeneities on
shortwave radiation: Comparison of cloud-top variability and internal heterogeneity. J. Atmos. Sci., 56, 4206–4224.
Villefranque, N., F. Couvreux, R. Fournier, S. Blanco, C. Cornet, V.
Eymet, V. Forest and J.-M. Tregan, 2019: A chain of Monte Carlotools to study interactions between radiation and clouds illustrated
on three applications. Submitted to J. Adv. Modeling Earth Sys.
Welch, R. M., and B. A. Wielicki, 1984: Stratocumulus cloud field reflected fluxes: the effect of cloud shape. J. Atmos. Sci., 41, 3085–
3103.
Wissmeier, U., R. Buras and B. Mayer, 2013: paNTICA: A fast 3D radiative transfer scheme to calculate surface solar irradiance for NWP
and LES models. J. Appl. Meteor. Climatol., 52, 1698–1715.
Wood, R., and P. R. Field, 2011: The distribution of cloud horizontal
sizes. J. Clim., 24, 4800–4816.
Yi, B., P. Yang, B. A. Baum, T. L’Ecuyer, L. Oreopoulos, E. J. Mlawer,
A. J. Heymsfield and K.-K. Liou, 2013: Influence of ice particle
surface roughening on the global cloud radiative effect. J. Atmos.
Sci., 70, 2794–2807.