Guided exercise: Energy Transport

In this guided exercise we will revisit the energy balance, this time with ERA5 data instead of CERES. The main novelty of this lesson is the computation of the meridional energy transport from the TOA, using the numerical methods we learned last week. The code is on me, the questions are on you! Like always, I recommend to download this notebook and run it by yourself instead of reading it online.

# Import the tools we are going to need today:
import matplotlib.pyplot as plt  # plotting library
import numpy as np  # numerical library
import xarray as xr  # netCDF library
import cartopy  # Map projections libary
import cartopy.crs as ccrs  # Projections list
# Some defaults:
plt.rcParams['figure.figsize'] = (12, 5)  # Default plot size

Data

We are going to work with a new file: ERA5_LowRes_AnnualAvg_eb.nc

ds = xr.open_dataset('../data/ERA5_LowRes_AnnualAvg_eb.nc')

Q: explore the data file. Some variables of the surface and TOA energy balance are missing, but they can be retrieved from the others: which ones?

ds

<xarray.Dataset> Size: 8MB
Dimensions:    (longitude: 480, latitude: 241)
Coordinates:
  * longitude  (longitude) float32 2kB -179.6 -178.9 -178.1 ... 178.9 179.6
  * latitude   (latitude) float32 964B 90.0 89.25 88.5 ... -88.5 -89.25 -90.0
Data variables:
    slhf       (latitude, longitude) float64 925kB ...
    sshf       (latitude, longitude) float64 925kB ...
    ssr        (latitude, longitude) float64 925kB ...
    str        (latitude, longitude) float64 925kB ...
    ssrd       (latitude, longitude) float64 925kB ...
    strd       (latitude, longitude) float64 925kB ...
    tisr       (latitude, longitude) float64 925kB ...
    tsr        (latitude, longitude) float64 925kB ...
    ttr        (latitude, longitude) float64 925kB ...

xarray.Dataset

• Dimensions:
  • longitude: 480
  • latitude: 241
• Coordinates: (2)
  • longitude
    (longitude)
    float32
    -179.6 -178.9 ... 178.9 179.6

    units :

    W m**-2

    long_name :

    longitude

    array([-179.625, -178.875, -178.125, ...,  178.125,  178.875,  179.625],
          dtype=float32)

  • latitude
    (latitude)
    float32
    90.0 89.25 88.5 ... -89.25 -90.0

    units :

    W m**-2

    long_name :

    latitude

    array([ 90.  ,  89.25,  88.5 , ..., -88.5 , -89.25, -90.  ], dtype=float32)

• Data variables: (9)
  • slhf
    (latitude, longitude)
    float64
    ...

    units :

    W m**-2

    long_name :

    Surface latent heat flux

    standard_name :

    surface_upward_latent_heat_flux

    [115680 values with dtype=float64]

  • sshf
    (latitude, longitude)
    float64
    ...

    units :

    W m**-2

    long_name :

    Surface sensible heat flux

    standard_name :

    surface_upward_sensible_heat_flux

    [115680 values with dtype=float64]

  • ssr
    (latitude, longitude)
    float64
    ...

    units :

    W m**-2

    long_name :

    Surface net solar radiation

    standard_name :

    surface_net_downward_shortwave_flux

    [115680 values with dtype=float64]

  • str
    (latitude, longitude)
    float64
    ...

    units :

    W m**-2

    long_name :

    Surface net thermal radiation

    standard_name :

    surface_net_upward_longwave_flux

    [115680 values with dtype=float64]

  • ssrd
    (latitude, longitude)
    float64
    ...

    units :

    W m**-2

    long_name :

    Surface solar radiation downwards

    standard_name :

    surface_downwelling_shortwave_flux_in_air

    [115680 values with dtype=float64]

  • strd
    (latitude, longitude)
    float64
    ...

    units :

    W m**-2

    long_name :

    Surface thermal radiation downwards

    [115680 values with dtype=float64]

  • tisr
    (latitude, longitude)
    float64
    ...

    units :

    W m**-2

    long_name :

    TOA incident solar radiation

    [115680 values with dtype=float64]

  • tsr
    (latitude, longitude)
    float64
    ...

    units :

    W m**-2

    long_name :

    Top net solar radiation

    standard_name :

    toa_net_upward_shortwave_flux

    [115680 values with dtype=float64]

  • ttr
    (latitude, longitude)
    float64
    ...

    units :

    W m**-2

    long_name :

    Top net thermal radiation

    standard_name :

    toa_outgoing_longwave_flux

    [115680 values with dtype=float64]

• Indexes: (2)
  • longitude
    PandasIndex

    PandasIndex(Index([-179.625, -178.875, -178.125, -177.375, -176.625, -175.875, -175.125,
           -174.375, -173.625, -172.875,
           ...
            172.875,  173.625,  174.375,  175.125,  175.875,  176.625,  177.375,
            178.125,  178.875,  179.625],
          dtype='float32', name='longitude', length=480))

  • latitude
    PandasIndex

    PandasIndex(Index([  90.0,  89.25,   88.5,  87.75,   87.0,  86.25,   85.5,  84.75,   84.0,
            83.25,
           ...
           -83.25,  -84.0, -84.75,  -85.5, -86.25,  -87.0, -87.75,  -88.5, -89.25,
            -90.0],
          dtype='float32', name='latitude', length=241))

• Attributes: (0)

Energy balance at the TOA

We’ll now compute the energy balance \(E_{TOA}\) at the top of the atmosphere (TOA), as well as its global average:

# constants
dx_rad = np.deg2rad(0.75)
coslat = np.cos(np.deg2rad(ds.latitude))
w = coslat / np.sum(coslat)
r = 6371000

eb_toa = ds.tsr + ds.ttr

ax = plt.axes(projection=ccrs.Robinson())
eb_toa.plot(ax=ax, transform=ccrs.PlateCarree())
plt.title('$E_{TOA}$')
ax.coastlines();

print('Incoming solar energy TOA: {:.0f} W m-2'.format((ds.tisr.mean(dim='longitude') * w).sum().data))
print('Average E TOA: {:.2f} W m-2'.format((eb_toa.mean(dim='longitude') * w).sum().data))

Incoming solar energy TOA: 340 W m-2
Average E TOA: 0.44 W m-2

Q: Is the energy budget for the total climate system balanced? Should it be? Compare the ERA5 radiative balance to the one we obtained with CERES data (0.97 W m\(^{-2}\))

Note: both datasets are in accordance with the current range of estimates of the Earth Energy Imbalance which ranges between 0.4 and 1 W m\(^{-2}\) (s1, s2)

Energy transport from the TOA Energy Balance

From the lecture, we remember that the total energy transport implied by the energy imbalance between the Equator and the Poles can be computed as:

\[ ET(\varphi)=\int\limits^{\varphi}_{-\pi/2}\int\limits^{2\pi}_{0} E_{TOA} R^2\cos\varphi'd\lambda d\varphi' \]

where \(R\) is the Earth radius (about 6371 km), \(\varphi\) the latitude and \(\lambda\) the longitude.

We are going to start by performing the inner integration:

\[\int\limits^{2\pi}_{0} E_{TOA} R^2\cos\varphi'd\lambda \]

Q: before I’ll do it for you, think about what we are going to compute: what is the purpose of \(\cos\varphi'\)? Of \(R^2\)?

intermediate_eb = (eb_toa * dx_rad).sum(dim='longitude') * r**2 * coslat

Q: what will be the shape of intermediate_eb?

We will now perform the outer integration:

\[ ET(\varphi)=\int\limits^{\varphi}_{-\pi/2} E (\varphi') d\varphi' \]

Where \(E (\varphi')\) is the intermediate variable we just computed.

# ERA5 goes from North to South but we integrate from S to N: invert
# I use cumsum here instead of trapz since the results are almost identical
et_toa = np.cumsum(intermediate_eb[::-1] * dx_rad, axis=0) / 1e15

et_toa.plot();
plt.hlines(0, -90, 90, linestyles=':')
plt.xlim([-90, 90]);

Q: Note that the transport is not zero at the North Pole. Can you guess why? Is this bad?

Energy balance at the surface

Now let’s compute the energy balance at the surface \(E_{Surf}\) as well as its average. Unlike with CERES data, we now have the opportunity to compute all fluxes at the surface (not only the radiative ones):

eb_surf = ds.str + ds.ssr + ds.sshf + ds.slhf
ax = plt.axes(projection=ccrs.Robinson())
eb_surf.plot(ax=ax, transform=ccrs.PlateCarree())
plt.title('$E_{Surf}$')
ax.coastlines();

print('Average SW Down: {:.0f} W m-2'.format((ds.ssrd.mean(dim='longitude') * w).sum().data))
print('Average SW UP: {:.0f} W m-2'.format(((ds.ssr - ds.ssrd).mean(dim='longitude') * w).sum().data))
print('Average LW Down: {:.0f} W m-2'.format((ds.strd.mean(dim='longitude') * w).sum().data))
print('Average LW UP: {:.0f} W m-2'.format(((ds.str - ds.strd).mean(dim='longitude') * w).sum().data))
print('Average Sensible HF: {:.0f} W m-2'.format((ds.sshf.mean(dim='longitude') * w).sum().data))
print('Average Latent HF: {:.0f} W m-2'.format((ds.slhf.mean(dim='longitude') * w).sum().data))
print('Average EB Surface: {:.0f} W m-2'.format((eb_surf.mean(dim='longitude') * w).sum().data))

Average SW Down: 188 W m-2
Average SW UP: -24 W m-2
Average LW Down: 339 W m-2
Average LW UP: -397 W m-2
Average Sensible HF: -17 W m-2
Average Latent HF: -83 W m-2
Average EB Surface: 6 W m-2

Q: Is the surface energy budget balanced? Should it be? Compare the values above with the values we discussed in the lecture. For reference, here is the figure again:

Note: errors in reanalysis data are aknowledged and discussed. See for example this paper for a discussion with ERA-Interim, the dataset prior to ERA5.

Energy balance of the atmosphere

I’ll now compute the energy balance of the atmosphere \(E_{Atm}\) as the residual of the energy balance at the surface and at the top of the atmosphere (note: this assumes that the entire residual of the surface EB is also taken from the atmosphere):

eb_atm = eb_toa - eb_surf
ax = plt.axes(projection=ccrs.Robinson())
eb_atm.plot(ax=ax, transform=ccrs.PlateCarree())
plt.title('$E_{Atm}$')
ax.coastlines();

print('Average E Atm: {:.2f} W m-2'.format((eb_atm.mean(dim='longitude') * w).sum().data))

Average E Atm: -5.73 W m-2

Q: Is the atmosphere energy budget balanced? Should it be?

The energy transported by the atmosphere can be computed with the exact same methodology:

\[ ET_{Atm}(\varphi)=\int\limits^{\varphi}_{-\pi/2}\int\limits^{2\pi}_{0} E_{Atm} R^2\cos\varphi'd\lambda d\varphi' \]

Let’s compute the energy transport the same way as we did before:

intermediate_eb = (eb_atm * dx_rad).sum(dim='longitude') * r**2 * coslat
et_atm = np.cumsum(intermediate_eb[::-1] * dx_rad, axis=0) / 1e15

et_toa.plot(label='$ET_{TOA}$')
et_atm.plot(label='$ET_{Atm}$')
plt.hlines(0, -90, 90, linestyles=':')
plt.xlim([-90, 90])
plt.legend();

Q Are you satisfied with the results? What are the major issues with this plot?

Energy transport of the atmosphere and the oceans

To compute the theoretical energy transport, we have to aknowledge that the energy budget has to be closed for the integration to work.

We now apply a bias correction to \(E_{TOA}\) and \(E_{Atm}\) by substracting the residual of the energy balance closure (so that the budget is balanced):

# Bias corrected TOA Transport
eb_toa_cor = eb_toa - (eb_toa.mean(dim='longitude') * w).sum()
intermediate_eb = (eb_toa_cor * dx_rad).sum(dim='longitude') * r**2 * coslat
et_toa = np.cumsum(intermediate_eb[::-1] * dx_rad, axis=0) / 1e15
# Bias corrected Atm Transport
eb_atm_cor = eb_atm - (eb_atm.mean(dim='longitude') * w).sum()
intermediate_eb = (eb_atm_cor * dx_rad).sum(dim='longitude') * r**2 * coslat
et_atm = np.cumsum(intermediate_eb[::-1] * dx_rad, axis=0) / 1e15

By noting that the total energy transport must be the sum of the transport by the atmosphere and by the oceans we can also compute \(ET_{Ocean}\):

et_ocean = et_toa - et_atm

et_toa.plot(label='$ET_{TOA}$', color='k')
et_atm.plot(label='$ET_{Atm}$', color='r')
et_ocean.plot(label='$ET_{Ocean}$', color='b')
plt.hlines(0, -90, 90, linestyles=':', color='k')
plt.xlim([-90, 90])
plt.legend();

Q: compare it to the plot by Fasullo and Trenberth (2008) that we analysed during the lecture.

Add-on: same with CERES data

ds = xr.open_dataset('../data/CERES_EBAF-TOA_Ed4.1_Clim-2005-2015.nc')
# constants
dx_rad = np.deg2rad(1)
coslat = np.cos(np.deg2rad(ds.lat))
w = coslat / np.sum(coslat)
r = 6371000
# Compute
eb_toa = (ds.solar_clim - ds.toa_sw_all_clim - ds.toa_lw_all_clim).mean(dim='month')
# Correct
eb_toa = eb_toa - (eb_toa.mean(dim='lon') * w).sum()
intermediate_eb = (eb_toa * dx_rad).sum(dim='lon') * r**2 * coslat
et_toa_ceres = np.cumsum(intermediate_eb * dx_rad, axis=0) / 1e15

et_toa.plot(label='$ET_{TOA}$', color='k')
et_toa_ceres.plot(label='$ET_{TOA_{CERES}}$', color='k', linestyle='--')
et_atm.plot(label='$ET_{Atm}$', color='r')
et_ocean.plot(label='$ET_{Ocean}$', color='b')
plt.hlines(0, -90, 90, linestyles=':')
plt.xlim([-90, 90])
plt.legend();

To go further

If you are interested, I recommend you have a look at the various terms of the TOA and surface energy balances and see if you can understand most of the patterns you observe. You can also plot their zonal average and analyse their respective importance according to latitude.