Top-of-atmosphere incident solar radiation for operational GraphCast

TL;DR. The library default for TOA solar radiation was a different physical quantity than GraphCast's training data. How I diagnosed and reconstructed the correct one.

Late 2023 I was wiring up the operational deployment of DeepMind’s GraphCast on KMA’s atmospheric model. The model expects a small set of input variables. One of them, toa_incident_solar_radiation (top-of-atmosphere incident solar radiation, accumulated J/m² over the previous hour), exists in the ERA5 reanalysis that GraphCast was trained against. KMA’s operational model (KIM) does not produce it directly. DeepMind’s guideline for running on non-ERA5 sources is to compute the variable analytically using the PySolar library. I did that and the result was noticeably wrong against ERA5. This post is about why, and how to fix it.

The variable, and why it has no operational source

toa_incident_solar_radiation is the solar flux arriving at the top of Earth’s atmosphere, integrated over a one-hour window, measured in J/m². It is a forcing variable in GraphCast: a deterministic function of (latitude, longitude, time) given a fixed model of Earth’s orbit. Because it does not depend on the atmospheric state, an operational NWP system (HRES, IFS, KIM) has no reason to forecast it. ERA5 carries it because ERA5 is reanalysis, which includes all forcing fields.

KIM does provide dswrtoa (downward shortwave radiation at TOA), but as an instantaneous value in W/m². Converting to J/m² needs an integral. A first attempt (multiplying the instantaneous value by 3600 seconds to get W·s/m² over the hour) ignores that the sun moves during that hour, and is biased everywhere there is a sun-edge.

DeepMind’s guideline, and the first failure

The official guidance is to compute toa_incident_solar_radiation via PySolar. PySolar exposes get_radiation_direct(when, altitude_deg) which returns the direct beam irradiance at the given sun altitude, in W/m². I wrapped this into hourly integration, tiled it across the (lat, lon) grid, and compared to ERA5.

Side-by-side with ERA5, the result was visibly off:

At the terminator (the boundary between day and night), the PySolar output was lower than ERA5 by close to an order of magnitude.
In the well-illuminated interior of the daylight hemisphere, PySolar was lower by ~15%.

Per-pixel ratio of the same fields. ERA5 / GDAPS-KIM (left) sits near 1 with the longitude offset visible at the sun-edge. The two ratios against PySolar's default (middle, right) show the terminator ring (deep red, ratio ≳ 10) and the ~15% interior under-prediction (light blue).

Both signs of error pointed at a consistent under-prediction: a physical term was missing.

PySolar’s `get_radiation_direct` returns the wrong physical quantity

Reading PySolar’s source revealed the issue. The library’s references describe two different quantities. First, the extraterrestrial solar insolation:

$I_0 = SC \cdot \left[1 + 0.034 \cos\left(\frac{360 n}{365}\right)\right] \quad (\mathrm{W/m}^2)$

where $SC$ is the solar constant (commonly 1377 W/m²) and $n$ is the day-of-year. Second, an attenuation model for atmospheric absorption (dust, water vapor, clouds, turbidity), treated as exponential decay along the slant path:

$I_B = A \, e^{-km}, \quad m = \frac{1}{\sin \beta}$

where $I_B$ is the beam portion reaching the surface, $A$ is “apparent” extraterrestrial flux, $k$ is the optical depth, and $m$ is the air-mass ratio for a sun altitude $\beta$ .

PySolar’s get_radiation_direct returns $I_B$ : the radiation that reaches the surface after atmospheric attenuation. ERA5’s toa_incident_solar_radiation is top of atmosphere, before any attenuation. The two differ by a multiplicative factor of $e^{-km}$ , which collapses toward zero where the air mass $m = 1/\sin \beta$ blows up: exactly the terminator region where I was seeing the worst error. The interior under-prediction is the same $e^{-km}$ factor at less extreme air mass.

The two are different physical quantities. PySolar’s default is correct for what a solar panel on the ground sees. It is the wrong quantity for toa_incident_solar_radiation.

Computing TOA from first principles

The correct TOA quantity is what hits a horizontal patch at the top of the atmosphere with no attenuation:

$I_{\text{TOA}}(\text{lat}, \text{lon}, t) = I_0(t) \cdot \sin\bigl(\beta(\text{lat}, \text{lon}, t)\bigr)$

when the sun is above the horizon, and zero otherwise. $I_0(t)$ is the extraterrestrial flux from the formula above; $\beta$ is the sun altitude at the given location and time, which PySolar exposes through its altitude utilities. The function is a few lines:

def get_tisr_map_v2(self, when):
    when = self.get_utc_when(when)
    altitude_map = self.get_altitude_map(when)        # degrees
    is_daytime = (altitude_map > 0)
    day = numeric.tm_yday(when)
    flux = 1377 * (1 + 0.034 * np.cos(2 * np.pi * day / 365))
    tisr = np.where(is_daytime,
                    flux * np.sin(np.radians(altitude_map)),
                    0.0)
    return tisr

This is a snapshot of TOA solar radiation at one instant, in W/m².

Hourly accumulation

ERA5’s toa_incident_solar_radiation is not a snapshot; it is an accumulation over the past hour, in J/m². To match, instantaneous snapshots have to be integrated:

$\text{TISR}(\text{lat}, \text{lon}, t) = \int_{t-1\text{h}}^{t} I_{\text{TOA}}(\text{lat}, \text{lon}, s) \, ds.$

A Riemann sum at 10-minute granularity is more than enough; the sun altitude does not change fast enough within an hour to need finer resolution:

def get_tisr_Jm2(self, when, step=60):
    """Accumulate W/m^2 over the past hour into J/m^2."""
    when = self.get_utc_when(when)
    start_when = when - datetime.timedelta(minutes=60)
    sub = datetime.timedelta(minutes=60 / step)
    dates = [start_when] + [start_when + sub * i for i in range(1, step)] + [when]
    tisr = np.zeros(self.longitudes_as_map.shape, dtype=np.float32)
    for date in dates:
        tisr += self.get_tisr_map_v2(date)
    tisr *= 3600.0 / step
    return tisr

The trailing 3600.0 / step converts “sum of W/m² snapshots” to J/m² over an hour.

Verification

The ratio plot ERA5 / (corrected PySolar) across the full lat-lon grid sat at ~1.0 everywhere, with no terminator artifacts. The same ratio against get_radiation_direct had shown the predicted pattern: a factor-of-ten under-prediction at the terminator and ~15% under-prediction in the daylight interior. The new method computes what ERA5 actually carries.

Ratio of TISR between sources. Corrected PySolar (left column) matches ERA5 and GDAPS-KIM at ~1 globally. PySolar's default get_radiation_direct (right column) under-predicts by an order of magnitude at the terminator and by ~15% in the daylight interior, against both reference fields.

The same custom integrator also resolved the KIM dswrtoa mismatch: ×3600 alone left a systematic 15°-longitude shift (the difference between the instantaneous KIM time stamp and the hour-accumulated ERA5 stamp), which the explicit hourly integration absorbs.

This fix shipped into the operational deployment in early 2024, predating DeepMind’s own patch addressing the same mismatch by a few months.

Closing

The lesson here is small and recurs in scientific ML wiring. Library defaults do not always describe the physical quantity they advertise; “solar radiation” routinely conflates what the sun emits and what reaches the ground, two quantities that differ by the air in between. The cost of a five-minute analytic check against ERA5 was the difference between a silent factor-of-ten bias at the terminator and an operational input that matched the training data.