Significance Early climate and weather models, constrained by computing resources, made numerical approximations on modeling the real world. One process, the radiative transfer of sunlight through the atmosphere, has always been a costly component. As computational ability expanded, these models added resolution, processes, and numerical methods to reduce errors and become the Earth system models that we use today. While many of the original approximations have since been improved, one—that the Earth’s surface and atmosphere are locally flat—remains in current models. Correcting from flat to spherical atmospheres leads to regionally differential solar heating at rates comparable to the climate forcing by greenhouse gases and aerosols. In addition, spherical atmospheres change how we evaluate the aerosol direct radiative forcing.

Abstract Sunlight drives the Earth’s weather, climate, chemistry, and biosphere. Recent efforts to improve solar heating codes in climate models focused on more accurate treatment of the absorption spectrum or fractional clouds. A mostly forgotten assumption in climate models is that of a flat Earth atmosphere. Spherical atmospheres intercept 2.5 W⋅m−2 more sunlight and heat the climate by an additional 1.5 W⋅m−2 globally. Such a systematic shift, being comparable to the radiative forcing change from preindustrial to present, is likely to produce a discernible climate shift that would alter a model’s skill in simulating current climate. Regional heating errors, particularly at high latitudes, are several times larger. Unlike flat atmospheres, constituents in a spherical atmosphere, such as clouds and aerosols, alter the total amount of energy received by the Earth. To calculate the net cooling of aerosols in a spherical framework, one must count the increases in both incident and reflected sunlight, thus reducing the aerosol effect by 10 to 14% relative to using just the increase in reflected. Simple fixes to the current flat Earth climate models can correct much of this oversight, although some inconsistencies will remain.

Our image of the atmosphere from space is that of a thin layer wrapped around the Earth: we see it as curved but thin enough so that locally, it can be treated as flat. In most all climate models, the 3D grid composing the atmosphere is tied to a 2D flat surface grid whose shape and area extends through the atmospheric layers to a subjective top-of-atmosphere height, typically 80 km for climate modeling. A simple picture of the spherical atmosphere and its flattened form is shown in Fig. 1. These 2 atmospheres are fundamentally different in how they intercept sunlight, and we identify 2 distinct factors: the light path and the spherical geometric expansion.

Fig. 1. Height (km) by distance (km) plot of surface (thick black line) and atmospheric layers (blue, every 5 km in height). (A) Spherical Earth and (B) flat Earth grid as used by most climate or Earth system models. The solid red lines show straight, unrefracted ray traces back to the sun for a SZA of 88°, which are (A) straight lines and (B) curved lines in the flat Earth grid. The dashed orange lines in B show the assumed ray trace to the sun in flat Earth models. For perspective on the curvature of the Earth, thin black lines demark 1° grid cells in the atmosphere: in A these are radial lines, and in B they are parallel lines. While not readily apparent on this scale, the length of the arc is 111 km at the surface and 112 km at 50 km height. We adopt the authalic radius of 6,371.0 km to ensure the same surface area as that of the true ellipsoidal Earth.

The light path is the most obvious error in the flat atmosphere and centers on the distortion of a straight line under the spherical-to-flat coordinate transformation. In Fig. 1A, the ray-tracing back to the sun at a surface solar zenith angle (SZA) of 88° is shown by the straight red lines originating at each layer. Here we assume parallel solar rays, ignoring sphericity associated with the sun–Earth distance (corrections are an order of magnitude smaller). The true light path ray traces in a flat atmosphere (Fig. 1B) appear as curves, reducing the zenith angle through each successive layer. The consequence of this is that the air mass factor, the ratio of path length to layer thickness, becomes smaller with altitude. In the flat atmosphere model where the air mass factors are constant, there is an exaggerated absorption of sunlight in the upper atmosphere, heating the upper atmosphere disproportionately relative to the lower.

Another aspect of the spherical light path is twilight. After the sun sets, we see the atmosphere overhead as sunlit. For example, when the sun is 4° below the horizon at the surface, the stratosphere is fully sunlit in clear skies with solar heating and active photochemistry, including Antarctic ozone depletion (1, 2). Thus, most atmospheric chemistry models use spherical solar ray tracing (3, 4), while most climate and Earth system models have retained the flat Earth ray tracing where atmospheric heating stops at an SZA of 90°. Several climate models use a curvature correction to SZA that limits the air mass factor (1/cos(SZA) in a flat atmosphere) approaching sunset but still treat the atmosphere as flat (Fig. 1B) and shut off the sun for SZA ≥ 90° (5⇓–7). The light paths in a spherical atmosphere can be simplified as straight lines (3) or can bend to include refraction (8, 9) (Methods). Flat atmospheres can include refraction, but at most this is invoked as a maximum air mass factor as SZA approaches 90°, which is still held constant through altitude.

In flat Earth models a column atmosphere is a polygonal cylinder with the same area A 0 from the surface to the top of atmosphere, while in a spherical atmosphere the columnar area expands like a cone as A 0 (1 + Z/R)2, with Z being the geometric height above the surface at radius R. Thus, a flat climate model is missing part of the climate system. The flat model solves the 1D planar hydrostatic equation with the gravitational acceleration fixed at its surface value g 0 ∼9.81 m⋅s−2. The mass of a layer in a flat column atmosphere with surface area A 0 is dM flat = A 0 dP / g 0 = – ρ ( P ) A 0 dZ pot , [1]where P is pressure, ρ is density, and Z pot is the geopotential height used in flat atmosphere models. Assuming a known profile of T(P) and composition as a function of P, then the profile of ρ(P, T) is also known. The profile Z(P) is integrated with the hydrostatic equation (right 2 terms in Eq. 1), again assuming that the gravitational acceleration g(Z) is a constant g 0 . The thickness of a layer in the flat geopotential atmosphere is less than that in the true geometric spherical atmosphere because g does not decrease with height, dZ pot = dZ / ( 1 + Z / R ) 2 . [2]This geometric–geopotential effect is noted in the National Center for Atmospheric Research (NCAR) Whole Atmosphere Community Climate Model (WACCM), which uses fixed g 0 , but it is dismissed as a problem for analyzing model output (10).

The spherical hydrostatic equation for a radially symmetric atmospheric layer with the same dP and ρ includes both conic expansion of the column and the correct gravitational acceleration to give a layer mass of dM sphere = ( 1 + Z / R ) 2 A 0 dP / g ( Z ) = ( 1 + Z / R ) 4 A 0 dP / g 0 = ( 1 + Z / R ) 4 dM flat . [3]The spherical mass is greater than the flat mass by a factor of 1 + 4H/R, where H ∼7 km is the density scale height (see equation 12 in ref. 11 and also ref. 12) This extra 1/2% of the atmosphere is optically active and intercepts more sunlight, bringing additional energy to the surface and climate system. The additional mass increases from 0% at the surface to 5% at a height of 80 km, with about half occurring in the 0 to 12 km height range. It cannot be corrected by simply decreasing g 0 .

A survey has found effectively flat atmospheres for the following climate models: Goddard Space Flight Center (Goddard Earth Observing System Chemistry Climate Model [GEOSCCM] and Modern-Era Retrospective analysis for Research and Applications, version 2 [MERRA-2]) (13, 14), Goddard Institute for Space Studies 2E (15, 16), Department of Energy (Energy Exascale Earth System Model [E3SM]) (17, 18), Geophysical Fluid Dynamics Laboratory (Atmospheric Model version 3) (19, 20), NCAR WACCM (10), and European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) (6, 21, 22). Recent development in terms of solar radiation has focused on more accurate treatment of water vapor (23), fractional clouds (21, 24), and averaging over time steps (6) but ignores the spherical nature of the atmosphere (12, 25).

Correcting the optical mass from flat to spherical presents a conundrum. We can calculate the absorption of solar energy in the climate system more accurately, but the energy deposited in each layer will be inconsistent with the mass of that layer in the flat atmosphere, resulting in excessive local heating. Unfortunately, the missing mass contributes to the heating of the climate system, including the surface. Here in Solar-J, we provide the capability for climate models to compute the scattering and absorption of sunlight in a spherical hydrostatic atmosphere that is consistent with their pressure and density profiles. It is up to the flat-atmosphere climate models to decide how to use this information.

When atmospheric layers expand with height, they increase the effective radius of the Earth in terms of capturing sunlight. For example, when calculating the sunlight incident on the climate system (i.e., stratosphere, troposphere, and surface; ozone, water vapor, aerosols, and clouds) with a refracting spherical expanded-mass atmosphere, we find about 2.5 W⋅m−2 more average radiation intercepted by the Earth than with a flat Earth model. If we include tropospheric aerosols and a large loading of stratospheric aerosols, the increase is about 3.2 W⋅m−2. This 0.7 to 1.0% increase is equivalent to extending the radius of the Earth by about 20 to 30 km. For example, even with clear-sky visible radiation, limb transmission of solar occultation measurements drops below 10% at 15 km height (26). This enhanced collection of sunlight is unevenly distributed geographically, however, with local increases of more than 10 W⋅m−2 in the twilight regions (7% of the globe at any given time). Photolysis rates that depend on overhead ozone column will increase at high SZA because of the reduced opacity of the spherical light path but will decrease at low SZA (high sun) because of the additional column mass using geometric heights. For example, at overhead sun, O 3 production from photolysis of O 2 in the lower stratosphere and upper tropical troposphere decreases by 4 to 6%, and excited atomic oxygen O(1D) in the troposphere decreases by 1%.

Solar-J version 7.6c, presented here, offers 4 options: option 0, flat, uses the flat atmosphere optical mass in each layer with a fixed air mass fraction of 1/cos(SZA) and a rigid cutoff at 90° SZA (Fig. 1B); option 1, sphr, adds true straight-line light paths (Fig. 1A) assuming the geopotential height grid, allowing for twilight; option 2, refr, allows these light paths to refract, extending surface sunset to SZAs of ∼90.5°; and option 3, geom, assumes increased optical mass in spherically expanded layers on a geometric height grid with refracted light paths. The solar intensity and flux deposited in each atmosphere are calculated accurately for each of these 4 models. For multiple scattering and absorption of diffuse solar radiation, Solar-J adopts the pseudospherical approach (8, 27) in which the solar intensity and flux are used in a 1D plane-parallel multistream Feautrier code (3, 28). The net solar flux deposited in each layer is derived from the divergence of the diffuse flux plus the direct beam. When invoking the spherical geometric atmosphere (geom), the direct beam’s deposition is enhanced by a factor (1 + Z/R)2 to account for the increase in atmospheric surface area. The increase in optical depth by the same factor due to the shift from geopotential to geometric heights has already been included in a corrected atmospheric profile of optical properties. These corrections for the 1D multiple-scattering ensure that all of the light scattered or absorbed within the expanding cone of atmospheric layers are included in that column’s energy budget. We do not consider the next order of correction for true spherical multiple scattering (29).

Profiles of the solar intensity for SZA from 0° to 94° are shown in Fig. 2A for the 4 models (flat, sphr, refr, and geom). This calculation assumes clear-sky conditions in the visible (600 nm) where the only optical components are O 3 Chappuis absorption and molecular Rayleigh scattering. The flat model starts deviating from sphr by SZA = 60°, and at SZA = 88°, flat overlaps with the SZA = 90° profiles of sphr, refr, and geom. The refracting codes (refr and geom) show increased intensity below 20 km compared with the straight-line light path sphr code for SZA > 90° as expected. The refr and geom differences, due only to the shift to geometric altitude, are barely discernible. The inflection point about 24 km in the SZA = 92° profile results from the solar rays at their tangent point passing below the weakly absorbing ozone layer. Similarly, when there are sharp cloud layers in the troposphere, there is often a full reversal with higher intensities in the lower layers for SZA > 90°. When the tangent point drops below the cloudy layer, the air mass factor for the cloudy layer is reduced (Fig. 1A). Such cases point to the difficulty in defining the incident solar flux in a spherical atmosphere.

Fig. 2. (A) Altitude profiles of the clear-sky intensity of sunlight at 600 nm as a function of SZA. Results are shown for flat grid (circles, stop at 88°), spherical light paths (sphr, dashed lines), spherical plus refraction (refr, solid lines), and spherical plus refraction plus geometric affects (geom, crosses). The inflection in the profiles for large SZA is caused by the ozone layer absorption. Below 30 km altitude, the flat grid model at 88° overlaps with the spherical models at 90°, consistent with spherical corrections that freeze the zenith angle at 88.6° until flat Earth sunset (6) (B) Relative change in incident solar flux at SZA = 88° caused by a 12-km liquid water cloud of OD = 4.2 (2 green colors) and by 0 to 5 km aerosol layer of OD = 0.08 (2 red-orange colors). Results are shown only for sphr (darker colors) and refr (lighter colors); effects of geom are small here. Change is relative to the clear-sky model for that light path assumption. The quantity plotted is change per kilometer so that the area under each curve is a measure of the total. The negative incident flux at the surface is plotted here (dashed lines) as if averaged over 20 km so that it can be compared with the solid curves. The summed net change in incident flux for atmosphere and surface is given. Optically thick clouds tend to have less incident flux than clear sky because of the great reduction of incident flux at the surface, while optically thin aerosol-like cloud tend to increase the incident flux. In all cases the total incident flux is greater than that of the flat grid model.

We define incident flux as that attenuated within the atmosphere or scattered or absorbed by the surface. For a flat atmosphere, it is simply and always equal to cos(SZA) times the solar flux. For a spherical atmosphere, it depends on the altitude and opacity of the atmosphere, and we show this in Fig. 2B for a case with 12-km-high cloud layer (optical depth = 4.2) and an aerosol-like layer in the lowest 5 km (optical depth = 0.08). Both cases use the same optical properties (liquid, stratus-like clouds) in an atmosphere with only weak Rayleigh scattering (800 nm). Both cloud layers (sphr = teal, refr = gray-green) and aerosol layers (sphr = red, refr = gold) intercept more sunlight, increasing the incident flux, but they also shadow the surface, decreasing the incident flux. The total atmosphere, surface, and net changes are shown in Fig. 2B. In this case and in general, optically thin aerosol or cloud layers have a net positive incident flux with the atmospheric capture exceeding the surface decrease, while optically thick layers reduce the incident flux. Nevertheless, even with optically thick clouds, the reduced incident flux remains still greater than that of the flat atmosphere.

Calculations shown here of the Earth’s solar radiation budget use the new Solar-J version 7.6c and meteorological data for January 2015 from the ECMWF OpenIFS system (cycle 38R1) (30, 31). We use 3-h averages for temperature, pressure, water vapor, ice-/liquid-water clouds on a 1.1° × 1.1° × 60-layer flat grid. Ozone is taken from a satellite climatology, and no aerosols are included in the base version. We integrate over 31 d with a 1-h time step. Increases in global incident solar radiation from successively implementing spher, refr, and geom are +1.6, +1.9, and +2.5 W⋅m−2, respectively, and are shown in Fig. 3 as a function of latitude. For perspective, a dashed line in Fig. 3 shows 1% of the total incident flux in the flat atmosphere. Errors in the full geom case are >2 W⋅m−2 over the sunlit globe and exceed 4 W⋅m−2 in 2 bands about 66°S and 66°N. For geom, about 60% of the increase in incident solar flux is absorbed by the climate system; this is split equally between atmosphere and surface, and in the atmosphere it is split equally between stratosphere and troposphere. Thus, the flat Earth climate models underestimate solar heating of the climate system by 1.5 W⋅m−2. See Table 1. For perspective, this difference is comparable to the total anthropogenic effective radiative forcing over the industrial era of 2.3 W⋅m−2 (32), which has produced a discernible climate change (33). Our geom correction is unlikely to effect the modeled climate change from preindustrial to present day, but it will likely shift both reference climates by amounts similar to the climate change, thus altering the model skill in reproducing the current climate (34).

Fig. 3. Increase in incident solar flux (W⋅m−2) relative to a flat Earth model as a function of latitude averaged over January 2015. Results are from a full solar heating code (24) updated to include the light path and geometric effects described here: sphr (blue), refr (red), and geom (green). See Table 1 notes, text, and SI Appendix. About 60% of the increase in incident solar flux is absorbed by the climate system, and 40% is reflected. The total incident flat Earth flux, scaled by a factor of 1/100, is plotted as a function of latitude (black dashed).

Table 1. Global mean solar fluxes (W⋅m−2) for January 2015 with clouds, no aerosols

Working in spherical atmospheres causes us to reexamine how we calculate the radiative effect of aerosols on climate. As one test case, we use the NCAR monthly mean 3D tropospheric anthropogenic aerosol climatology for surface area active in heterogeneous chemistry (35). When converted to a mix of aerosol types, we calculate a global mean 600-nm optical depth of 0.088, typical of others (36). For January we calculate the aerosol direct radiative forcing (also known as instantaneous radiative effect). Our aerosols increase reflected sunlight by 1.56 W⋅m−2, but they reduce the total heating by only 1.35 W⋅m−2 because they catch more sunlight, increasing the incident flux by 0.21 W⋅m−2 (Table 2). The balance between atmospheric heating and surface cooling here is similar to recent model comparisons (37), and we find that the flat atmosphere high-bias errors of about 7% are comparable to the identified model errors (figure 2 of ref. 37). In another test case, we take the stratospheric aerosol loading and effective radius from the geoengineering model intercomparison project (38, 39). Our 600-nm optical depth of about 0.094 is similar to other models in the project and also to the levels following large volcanic eruptions. For stratospheric sulfate aerosols, our test is more rigorous because we can accurately calculate and resolve their optical properties and effects through 1) the index of refraction for different weight percent solutions of sulfuric acid, 2) Mie scattering, and 3) 8-stream multiple scattering. The geoengineered cooling derived from the increase in reflected sunlight, −2.8 W⋅m−2, must be offset by the increase in incident flux, +0.5 W⋅m−2. Thus we need to evaluate aerosol–climate effects in terms of the change in heating of the climate system rather than by the change in reflected solar energy.

Table 2. Aerosol direct radiative forcing (W⋅m−2) including incident flux for January 2015 all-sky conditions

Conclusions The approach here is only an approximation to radiative transfer solutions for a spherical hydrostatic atmosphere, but it is a low-cost option to greatly improve the accuracy of photolysis, photosynthesis, and solar heating in current models. Calculation of the direct solar beam is rigorous. Reverting to a plane-parallel scattering code does not address the 3D nature of scattered light (e.g., sides of clouds (22) and shift in quadrature angles with altitude (Fig. 1)), but it is still an improvement over flat atmospheres. Implementing these spherical corrections within an existing flat solar heating code should be straightforward and without significant penalty on computational efficiency. In Solar-J, the ray tracing and geometric corrections are calculated once for each column atmosphere and SZA and then applied to all wavelengths and combinations of cloud cover. This is an approximation as refraction depends somewhat on wavelength and atmospheric scale height. The attenuation of solar flux uses the optical properties for each wavelength band and applies Beer’s Law along each ray path to the sun. The solar intensity at each layer edge has a unique ray path and is not directly related to the intensity at the layer above. A possibly disruptive effect of spherical atmospheres is that a radiation call must be made for about 58% of the grid cells that are sunlit instead of exactly 50% and thus may affect load balancing. Most chemistry models have pseudospherical atmospheres with twilight photochemistry and have dealt with load balancing across column atmospheres. Calculation of solar energy absorption in a spherical hydrostatic atmosphere like the Earth’s using the geom option would seem to be the correct approach for climate modeling as the corrections relative to a flat, reduced mass atmosphere penetrate through the atmosphere to the ocean and land systems. Nevertheless, it raises questions of consistency and how to distribute the energy within the atmosphere. Our results here are based on 1-h model stepping and produce the same magnitude wavenumber 24 oscillations in incident flux and atmospheric heating found previously (6, 40). Using a spherical atmosphere reduces the amplitude of both errors by about half because it softens the transition to dusk but does not eliminate the need for some approach to averaging the heating rates over the full time step at each longitude (6). In contrast with ref. 6, our code has no obvious bias in stratospheric heating rates because we do not alter the air mass factors approaching sunset: differences between 1-h and 30-min time steps is less than 0.001 K⋅d−1. The difference in stratospheric heating rates with sphericity is quite large, however: +0.6 K⋅d−1 near the stratopause in 2 bands about 66°S and 66°N. Spherical atmospheres should also affect the thermal infrared cooling of the Earth. Effects would include the increase in overall opacity (thicker geometric layers) but also smaller opacity along nonzenith ray paths to space (air mass factors decrease with height along a ray path; Fig. 1A). The sign of the bias is not obvious, but an investigation using the tools developed here could be readily made. For planets with atmospheric scale height to radius ratios greater than the Earth’s (∼0.1%), such as Venus or Mars (∼0.3%), spherical and geometric effects will become more pronounced. Without implementing these spherical corrections in a climate model, we can only speculate on how they would affect the simulated climate. Does an additional global 2 W⋅m−2 out of a total heating of 240 W⋅m−2 matter? Since climate models can clearly detect the climate change caused by similar levels of human forcing (41), the shift in climate caused by flat atmosphere errors should be obvious. Also, 2 W⋅m−2 is typical of the model-to-model differences in atmospheric heating (42). The spherical corrections, even for the flat atmosphere mass, are not uniform and thus may cause much larger corrections in the stratosphere and winter poles. The spherical photochemistry modules already being used in Earth system models can be merged with the solar heating modules to maintain physically consistent codes.

Methods The solar heating rates calculated here use updated versions of the standalone Cloud-J code (43) and the Solar-J code (24) embedded within the University of California, Irvine, chemistry-transport model. This latter model was run for January 2015, taking hourly steps and using meteorological data from the ECMWF IFS forecasts at T159L60 resolution (about 1.1° × 1.1° resolution with 60 layers; for details, see ref. 31). Ozone profiles are taken from a satellite climatology included in Cloud-J and Solar-J. Solar-J uses the superbins and water vapor subbin absorption of Rapid Radiative Transfer Model for shortwave radiation (RRTMG-SW) for wavelengths greater than 600 nm. Cloud-J includes the same superbins and cloud absorption as Solar-J but does not include water vapor absorption. The new versions 7.6c of both models in standalone forms are provided in ref. 44. The stratospheric geoengineering sulfate aerosols are taken from The Geoengineering Model Intercomparison Project (GeoMIP) specification for experiment G4-SSA (38). Data are taken from the midscenario, near–steady-state values for zonal mean effective radius and H 2 SO 4 density in the file: geomip_ccmi_2020-2071_volc_v3.nc, from the GeoMIP website. The density was converted into a more general form of mass ratio (kg H 2 SO 4 per kg air) and scaled up to include the 25 wt % of H 2 O in the aerosols. The first 8 terms of the scattering phase function, the single scattering albedo, and the extinction efficiency (often denoted by Q) were calculated with a Mie code assuming a refractive index for 75 wt % sulfuric acid aerosols (density = 1.746 g⋅cm−3), which is typical for the stratospheric sulfate layer. The aerosol size distribution was assumed to be log-normal with a width of 1.25. Separate calculations for the 27 wavelength bins of Solar-J were made for effective radius of 0.02, 0.04, 0.08, 0.10, 0.14, 0.20, 0.30, 0.50, 0.75, 1.0, 1.2, 1.4, 2.0, 3.0, and 5.0 microns. For the GeoMIP aerosols, we selected the optical properties from the one of these 15 precalculated size bins with the effective radius closest to that prescribed in the layer. We calculate an annual mean sulfur loading of 2.08 TgS, consistent with the GeoMIP documentation (38). For January 2015, we have a mean optical depth at 600 nm of 0.0935. GeoMIP also causes photochemical perturbations. In the lower stratosphere, O 3 production from photolysis of O 2 decreases by 2 to 3% while atomic O increases by 1 to 4%, thus generating a direct photochemical depletion of ozone independent of heterogeneous chemistry. In the tropospheric boundary layer, O(1D) decreases by 1.3% reducing hydroxyl radical (OH) concentrations, but this does not include the expected decreases in stratospheric ozone. For this GeoMIP calculation, we believe our direct radiative results to be accurate, more so than most of the radiation codes used in the GeoMIP models: the optical properties are based directly on the refractive index for sulfuric acid–water mixtures and averaged over the wavelength bins, we include the wt % of H 2 O, and the multiple scattering solution uses 8-streams for the diffuse radiation with no equivalent–isotropic reduction in optical depth. The tropospheric aerosols are taken from 3D monthly mean, 10-y averages of surface area densities (µm2 per cm3) from the Community Earth System Model version 1.2 (35) supplied by S. Tilmes on 17 February 2016. The surface area includes the wet radius. Being intended for heterogeneous surface chemistry, it does not include sea salt, dust, or specific hydrophobic aerosols but is primarily biomass burning aerosols, industrial sulfate, and other hydrophilic aerosols. For optical properties, we split the surface area into 2 types: 1) combustion aerosol with effective radius of 0.14 µm, single scattering albedo at 600 nm = 0.89, and extinction efficiency at 600 nm Q = 1.06 and 2) dust aerosol with effective radius of 0.15 µm, single scattering albedo = 0.98, and extinction efficiency Q = 1.28. These 2 were chosen because they are available as standard aerosol types in Fast-J (#9 and #11, respectively) and they give typical values for overall single scattering albedo, Angstrom coefficient, and optical depth at 600 nm of 0.0877. We recognize that dust aerosol was not included as anthropogenic aerosol in the CESM data. In this paper, the tropospheric aerosol loading is chosen as a didactic example for anthropogenic aerosols, to demonstrate the role of spherical atmospheres. A more accurate evaluation of the direct radiative forcing of tropospheric aerosols would have to include the full range of anthropogenic aerosol optical types and sizes. The treatment of air mass factors in a spherical refracting atmosphere are described in the equations A.1–A.34 of SI Appendix. The radiative transfer solution is done for each grid cell column atmosphere, assuming that it is horizontally homogeneous in spherical shells. The unrefracted SZA for each column atmosphere is the same at each layer edge throughout the column and is determined by latitude, longitude, solar declination, and local time. With refraction, the zenith angle of the refracted beam is reduced relative to that of the unrefracted solar beam, in proportion to vertical gradients in density. We adopt 1.00030 as the index of refraction of air for all wavelengths at the surface with a fixed density scale height of 8 km. This assumption can be revised to include realistic density variations and wavelength dependence of the refractive index (e.g., from 1.00030 in the UV to 1.00027 in the near infrared). Presently, Solar-J is designed for the stratosphere and troposphere, and thus, we shut off the heating rate calculation for SZAs greater than 98°, corresponding to no direct solar beam below 62 km. The absorption and extinction of the direct solar rays are calculated with correct spherical geometry, using geopotential (sphr and refr) and geometric (geom) heights. This calculation provides the solar source terms for the scattering calculation, for which we revert to our plane-parallel (flat atmosphere), multiple-scattering code (3), increasing the optical mass of the flat atmosphere layer in the geom option to account for the additional thickness and areal extent of the layer. This method of combining spherically calculated direct solar terms with a 1D plane-parallel model for the diffuse solar radiation provides a reasonable approximation to spherical atmosphere that is exact in limiting cases such as optically thin layers. The standalone codes for Solar-J and Cloud-J versions 7.6c, along with other data in the paper are published with University of California, Irvine Dash (44).

Acknowledgments This work was developed with support from the US Department of Energy, Office of Science, Biological and Environmental Research Program (award DE-SC0012536); Lawrence Livermore National Laboratory (subcontract B628407) under the E3SM project; and the NASA Modeling, Analysis and Prediction program (award NNX13AL12G). We greatly appreciate the efforts of Tim Cronin and the anonymous reviewers for improving the overall presentation and identifying our inconsistent treatment of the spherical hydrostatic atmosphere, which is now corrected.

Footnotes Author contributions: M.J.P. designed research; M.J.P. and J.C.H. performed research; M.J.P. and J.C.H. analyzed data; and M.J.P. and J.C.H. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1908198116/-/DCSupplemental.