Studying climate of the Maastrichtian (~70 Million years before present)

Note

This document represents one of seven different tutorial handouts. Each handout is similar regarding the analyses that shall be performed on climate data. Differences in handouts relate to the climate states that are being investigated and to the climate model employed to create the climate data.

A very short introduction to climate models

In this tutorial we will study selected characteristics of a climate state that is significantly warmer than the one we are currently experiencing. We analyze data that has been produced by a numeric climate model (in the following termed climate model for brevity) and quantify various climate-related physical characteristics. Climate models are comparable to weather prediction models that are routinely used to create weather forecasts. Climate models employ a similar formalism and model structure as weather forecast models. An aspect in which both types of models differ is that their use is based on different methodologies: meteorologic conditions in the short term future (e.g. over the coming days) are mostly driven by the initial state of the atmosphere and by chaotic behavior that characterizes the trajectory of the state in the weather forecast model's phase space. In contrast, for the analysis of climate, that is defined as the long-term (~30 years, for example) time average of weather conditions, the initial state of the climate system is less important. When using climate models we rather focus on the impact of boundary conditions (e.g. strength of solar irradiation at the top of the atmosphere, composition of the atmosphere, distribution of land masses) and study their impact on the state of the climate system. For practical considerations a focus on boundary conditions is actually an advantage. In comparison to initial conditions, that are needed for weather forecast models (e.g. spatial distribution of temperature at 12am UTC), our knowledge of climate system boundary conditions is much more precise - both for past and present. This explains why it is useful to employ reconstructions of relevant characteristics of climate drivers to study a climate state that was present several million years ago whereas it is completely unrealistic to aim for correctly forecasting weather conditions in, lets say, four weeks. The latter is practically impossible due to uncertainties in the initial condition and chaotic behavior of weather. In climate research such problems are largely absent as short-term fluctuations in weather conditions are to a large degree averaged out.

Depending on the handout of the tutorial you have been provided with the climate data at hand has been produced employing either the COmmunity Earth System MOdels (COSMOS) or the Max Planck Institute Earth System Model (MPI-ESM). These modelling systems are capable of simulating the impact of changes in boundary conditions on the states of atmosphere, land surface, and ocean. In the context of this tutorial the most relevant boundary condition is the amount of carbon dioxide available in the atmosphere. Other boundary conditions that may be changed in simulations with these models (and that are changed in simulations of some of the tutorial handouts), are land surface characteristics, including height and extent of ice sheets, state of ocean gateways, location and elevation of mountain chains, etc. Details of model setup and application of a climate model for paleoclimate research are described by Stepanek et al. (2020) at the example of the COSMOS.

When interpreting results of climate simulations it is important to consider that not all components of the real Earth System are considered in the model. Per definition of a model also climate models will always provide imperfect representations of reality. When interpreting results of climate simulations it is hence of profound importance to know

• which components of the real Earth System are explicitly resolved by the model and
• in which respect characteristics of the Earth System are strongly simplified in the model for practical reasons.
  

In the model setup of the climate model, that has been employed for creating the climate simulations analysed here (refer to Figs. 1 or 2), circulation of atmosphere and ocean, and the most relevant climate-related processes including sea ice dynamics, are resolved. In addition, selected climate-related processes that occur at the land surface are considered in the model. Yet, there is a multitude of compartments of the Earth System not resolved by the model. This problem is in the following discussed at the example of ice sheets. Ice sheet dynamics are not dynamically resolved by the model setup employed here. This means that ice sheets cannot build up or decay during a simulation with the model as they would do in reality in response to changes in ambient climate. Rather, in the simulations that are analysed here, the state of an ice sheet is considered as being a part of the boundary conditions that drive the model. To this end, both height and extent of ice sheets are prescribed to the model based on auxiliary information that is provided, for example, by means of ice sheet simulations conducted for the relevant time period in past or future. There are two consequences of this approach: first, simulated climate characteristics will always consider the state of a fixed prescribed ice sheet; second, climate and land ice may not necessarily be in equilibrium with each other.

Fig. 1: Components of the Community Earth System Models. Shown are the climate-subsystems resolved by the model - these are atmosphere, ocean, and land surface. The schematic specifies names of component models as well as their horizontal and vertical resolution. The model component OASIS represents the coupler that facilitates that fluxes of energy, mass, and momentum are exchanged between atmosphere and ocean - without coupling such fluxes atmosphere and ocean could not react to changes in each other's states. A result would be that important patterns of climate variability, e.g. El Nino Southern Oscillation, that are caused by communication between atmosphere and ocean, would not occur in simulations performed with the model. Source: Stepanek et al. (2020).

Fig. 2: Components of the Max Planck Institute Earth System Model. Shown are the climate-subsystems resolved by the model - these are atmosphere, ocean, and land surface. The schematic specifies names of component models as well as their horizontal and vertical resolution. The model component OASIS represents the coupler that facilitates that fluxes of energy, mass, and momentum are exchanged between atmosphere and ocean - without coupling such fluxes atmosphere and ocean could not react to changes in each other's states. A result would be that important patterns of climate variability, e.g. El Nino Southern Oscillation, that are caused by communication between atmosphere and ocean, would not occur in simulations performed with the model. Source: Stepanek (unpublished).

Aims of the tutorial

The tutorial aims at providing to you the opportunity to perform analyses of climate model output as they are routinely done in climate research. One aim of the tutorial is to perform this exercise in a setting that is as realistic as possible from the viewpoint of climate research. To this end we employ in this tutorial standard research tools. The models that have been employed to create the climate data provided to you are actively used in climate research. The climate states, respectively climate epochs, considered in the tutorial have not been specifically selected for teaching purposes - rather, they are a subject of active research in paleoclimate. You have been provided with model output that would also be used by us for preparing scientific publications. Most of the climate data employed in the tutorial actually has been published in peer-reviewed journals.

In a first step of our tutorial we will prepare graphics that illustrate some characteristics of the analyzed climate states. These will encompass geographic maps of time-average quantities on the one hand and time series of spatial-average quantities on the other hand. Derived illustrations will enable us to gain a first idea of

• how climate differs between experiment and reference in various regions of the world and
• how large climate variability is over time.

In a second step we will derive first order climate characteristics. We will determine how sensitive the simulated global surface air temperature is to a change in carbon dioxide. This will be done at the example of a doubling of carbon dioxide from Pre-Industrial conditions. For climate states where boundary conditions in the model are nearly identical, the only exception being differences carbon dioxide concentrations, this quantity is termed Climate Sensitivity. For climate simulations where boundary conditions differ in more characteristics than carbon dioxide alone (including, for example, the extent and height of ice sheets) this quantity is rather termed Earth System Sensitivity. This reflects that in such simulations not only climate reacts to carbon dioxide, but the Earth System as a whole is changed, including for example height and extent of ice sheets. Furthermore, we will quantify the amplitude of Polar Amplification. This term describes the tendency of the climate system to amplify temperature change at high latitudes. A modest low-latitude or global average warming may be expressed as substantial temperature change in polar regions. This effect is a common feature of warm climate states and also visible in the model data analysed in the handout to this tutorial (Figs. 4. and 5). By the end of the century the global average temperature change due to anthropogenic emissions will be exceeded by a factor of two to three in the Arctic (Davy and Outten, 2020), with far reaching consequences for a vulnerable cryosphere in polar regions. A bit of background on these two important climate metrics is provided below:

• Climate Sensitivity (Earth System Sensitivity; Lunt et al., 2012) quantifies the susceptibility of global average near surface air temperature to perturbations in carbon dioxide. Climate Sensitivity (Earth System Sensitivity) is defined as the global average temperature change that arises if carbon dioxide is doubled from the Pre-Industrial level (280 parts per million by volume, ppmv), i.e. to 560 ppmv. This quantity is hence of utmost interest if we would like to estimate climate warming that is to be expected as a result of humankind perturbing greenhouse gas concentrations in the atmosphere. The concepts of Earth System Sensitivity and Climate Sensitivity differ from each other in that they consider different types of feedbacks in the climate system. Climate Sensitivity only takes into account fast feedbacks in atmosphere and ocean. These include water vapour, sea ice, and clouds. Yet, on long time scales also slow feedbacks become relevant in shaping characteristics of a climate state. These include ice sheets, vegetation, and surface albedo who all play a role in driving climate on the long term. Such slow feedbacks are additionally considered by extending the concept of Climate Sensitivity to that of Earth System Sensitivity. Whether a set of simulations enables quantification of Earth System Sensitivity depends on the boundary conditions that were used to force the climate model. In our analysis Earth System Sensitivity can be considered if we study climate states that have been simulated employing model configurations characterized by differing Earth surface characteristics, whereas Climate Sensitivity is quantified if the land surface conditions are identical. For example, Climate Sensitivity can be quantified by comparing temperatures derived from a set of two simulations that were both forced with modern land surface conditions, where one simulation is forced with the Pre-Industrial concentration of carbon dioxide and the other simulation is forced with a concentration of 560 ppmv of carbon dioxide. In contrast, Earth System Sensitivity may be studied if the simulation with 560 ppmv of carbon dioxide also features land surface characteristics indicative of respectively elevated surface temperatures - for example, a representation of ice sheet extent and height that is representative for a cryosphere under the impact of carbon dioxide concentrations doubled from Pre-Industrial conditions.
• Polar Amplification (e.g. Pithan and Mauritsen, 2014) describes the fact that high latitudes of the Earth System generally warm much more than low latitudes. Polar Amplification is defined as the ratio of warming in the region north of 60°N (south of 60°S) with respect to warming at low latitudes (0-30°) of the respective hemisphere. Polar Amplification is defined for both hemispheres and its value may differ between northern and southern polar regions. In this tutorial we focus on amplification of warming in the polar region of the Northern Hemisphere, hence on Arctic Amplification.

When analysing the climate state in this tutorial we will do so based on an anomaly approach. That means that we compare a climate simulation of interest to a reference state. This approach enables us to relate specific changes in the climate state between the simulation of interest and the reference state to respective changes in boundary conditions. In climate modelling the state of the Pre-Industrial (~1850 CE) is often taken as a reference. Rationale is that the Pre-Industrial provides a more stable climate state than our current ambient climate that has strong transient characteristics driven by ever increasing levels of greenhouse gases of anthropogenic origin.

Characteristics of the climate model output to be analyzed

Climate model output is available as a spatially resolved time series of monthly mean output. The time series are available over 100 model years. Their length is hence sufficient towards enabling the analysis of climate characteristics in contrast to properties of more short term weather. There are two climate simulations considered per tutorial handout. These climate states differ in the concentration of carbon dioxide and/or earth surface characteristics considered in the model. Comparing two simulations allows quantification of the impact of specific differences in boundary conditions on the simulated climate state. Data is available for the whole global domain so that it is possible to compute spatial averages over smaller regions. While both the COSMOS and the MPI-ESM simulate both atmosphere and ocean quantities, here we focus on the output of the atmosphere model as the atmosphere is the climate realm for which quantities like Climate Sensitivity, Earth System Sensitivity, and Polar Amplification are estimated. For COSMOS atmosphere model output has a spatial resolution of 3.75°x3.75° (longitude/latitude), for MPI-ESM the resolution is 1.88°x1.88° (longitude/latitude). In both cases spatial resolution of the model is rather coarse if compared to that of weather forecast models. Yet, for paleoclimate research the resolution is in the normal range. In climate research, and in particular in paleoclimate research, where slow model components like the ocean have to be brought into equilibrium with climate forcings that significantly differ from those of today, coarse model resolution of long-term climate simulations as compared to high resolution of short-term weather simulations is a necessity. Computational power is limited. Hence, if one would like to create model simulations over a longer time period then one has to accept that this normally will come at the cost of a lower spatial resolution of the model.

Fig. 3: Time line of the Cenozoic (since about 66 Million years ago). Top: Subdivision of the Geologic Era of the Cenozoic into Periods and Epochs following the time scale of the International Commission on Stratigraphy. Middle: Climate evolution across the Cenozoic depicted by anomalies of global average surface temperature at a specific time with respect to the global average surface temperature of the Pre-Industrial reference period. Bottom: Reconstructed past and projected future concentration of carbon dioxide in the atmosphere. Time line of geologic / climatic history from left to right. Note the use of different scales on the time axis. The rightmost section of the graph depicts projected evolution of global average temperatures until 2300 CE. Similarity of past and future global average temperatures illustrates that the climate of the Cenozoic is relevant for studying potential near future conditions. This inference is a motivation for conducting paleoclimate research. Source: ECORD/IODP (https://twitter.com/ECORD_IODP/status/1407018887523278854/photo/1).

The Maastrichtian - the dinosaurs "last hour"

The Maastrichtian (~72.1 - 66.0 Million years before present; Fig. 3) is a geologic age that directly precedes the Cenozoic, Earth's current geological era. The Maastrichtian is the last age of the Cretaceous, the end of which is marked by the Cretaceous-Paleogene extinction event. This event marks the extermination of a significant number of species of the animal kingdom, including a large number of dinosaurs. Extinction was triggered by the impact of an asteroid. Global climate was not as hot and humid anymore as it had been during earlier stages of the Cretaceous, yet in comparison to later parts of the Cenozoic global temperatures were still extreme. Noteworthy is that the Maastrichtian was still so warm that ice sheets were absent in the Earth System.

When studying the Maastrichtian we have to consider a geography that appears very unusual from the viewpoint of today (Müller et al., 2016; Niezgodzki et al., 2017). The distribution of continents was still very different from current geographic conditions. The Atlantic Ocean, that separates the Americas from Europe and Africa, already existed. Yet, during the Maastrichtian the ocean basin was still rather narrow, and a deep ocean section, separating Europe from North America, was still absent. Furthermore, the energy output of the sun was in comparison to today smaller, which resulted in a solar constant that was about 1% smaller than today (Niezgodzki et al., 2017, and references therein). From that point of view the Maastrichtian is a less suitable analogue to potential future climate conditions than other periods of the Cenozoic, like Pliocene or Eocene (Burke et al., 2018), where climate is warmer than at present but at similar geography and solar constant. Yet, the Maastrichtian still provides a useful laboratory for climate science, for example towards the study of the evolution of meridional circulation and related northward heat transport in the Atlantic Ocean basin. This type or circulation established after the opening of a connection between Atlantic and Arctic Ocean and lead to significant warming of large regions of Northern Hemisphere high latitudes (Hunter et al., 2013). Studying evolution and effect of meridional circulation in the Atlantic Ocean basins is of profound interest as the related heat transport is today of major importance for the expression of high latitude climate in the Atlantic/Arctic domain.

In this tutorial we study a Maastrichtian climate state that is driven by a reconstruction of land surface and ocean bathymetry conditions that are representative for that time period. In the Maastrichtian climate simulation the following boundary conditions are modified from those of the Pre-Industrial reference state: ocean gateways, ocean bathymetry, land surface elevation, ice sheet extent, as well as vegetation and related land surface characteristics. In addition there is a change of solar constant (1% reduction with respect to today). The employed carbon dioxide concentration is double that of the Pre-Industrial. The design of the simulation enables us to study a fundamental climate characteristic during the Maastrichtian: Earth System Sensitivity in a world that is bare of any ice sheets. Studying Earth System Sensitivity is possible in this setup as the simulation of the Maastrichtian does also consider changes of land surface characteristics that reflect a warmer climate state. These have been reconstructed from various sources of (mostly geologic) information.

For simplicity, we refer to the simulation with Maastrichtian boundary conditions and 560 ppmv of carbon dioxide as experiment. Similarly, we refer to the Pre-Industrial reference state, that is driven by a modern configuration of earth surface conditions and by a Pre-Industrial concentration of carbon dioxide (280 ppmv), as reference.

A detailed description of simulated climate states employed here is given by Niezgodzki et al. (2017). In case you would like to get additional information note that simulation experiment of this tutorial is in the publication by Niezgodzki et al. (2017) named C-560, whereas simulation reference of this tutorial is in their work named PI.

References

Burke, K. D., Williams, J. W., Chandler, M. A., Haywood, A. M., Lunt, D. J., and Otto-Bliesner, B. L.: Pliocene and Eocene provide best analogs for near-future climates, P. Natl. Acad. Sci. USA, 115, 13288–13293, https://doi.org/10.1073/pnas.1809600115, 2018.

Davy, R. and Outten, S.: The Arctic Surface Climate in CMIP6: Status and Developments since CMIP5, Journal of Climate, 33, 8047–8068, https://doi.org/10.1175/JCLI-D-19-0990.1, 2020.

Hunter, S. J., Haywood, A. M., Valdes, P. J., Francis, J. E., and Pound, M. J.: Modelling equable climates of the Late Cretaceous: Can new boundary conditions resolve data–model discrepancies?, Palaeogeogr. Palaeoclimatol. Palaeoecol., 392, 41–51, https://doi.org/10.1016/j.palaeo.2013.08.009, 2013.

Lunt, D., Haywood, A., Schmidt, G., Salzmann, U., Valdes, P. J., and Dowsett, H. J.: Earth system sensitivity inferred from Pliocene modelling and data, Nature Geosci., 3, 60–64, https://doi.org/10.1038/ngeo706, 2010.

Müller, R. D., Seton, M., Zahirovic, S., Williams, S. E., Matthews, K. J., Wright, N. M., Shepard, G. E., Maloney, K. T., Barnett-Moore, N., Hosseinpour, M., Bower, D. J., and Cannon, J.: Ocean basin evolution and global-scale plate reorganization events since Pangea breakup, Annu. Rev. Earth and Planet. Sci., 44, 107–138, https://doi.org/10.1146/annurev-earth-060115-012211, 2016.

Niezgodzki, I., Knorr, G., Lohmann, G., Tyszka, J., and Markwick, P. J.: Late Cretaceous climate simulations with different CO2 levels and subarctic gateway configurations: A model-data comparison, Paleoceanography, 32, 980–998, https://doi.org/10.1002/2016PA003055, 2017.

Pithan, F. and Mauritsen, T.: Arctic amplification dominated by temperature feedbacks in contemporary climate models, Nature Geosci., 7, 181–184, https://doi.org/10.1038/ngeo2071, 2014.

Stepanek, C., Samakinwa, E., Knorr, G., and Lohmann, G.: Contribution of the coupled atmosphere–ocean–sea ice–vegetation model COSMOS to the PlioMIP2, Clim. Past, 16, 2275–2323, https://doi.org/10.5194/cp-16-2275-2020, 2020.

Analyses

loading Python functionality necessary to perform the analysis

#provide auxiliary functionality by loading specific packages
from netCDF4 import Dataset #provides functionality to read/write binary NetCDF data
import matplotlib.pyplot as plt #provides functionality for graphical representation of data
import matplotlib.colors as colors #provides functionality for optimizing plot colors
import numpy as np #provides functionality for scientific computing
import cartopy.crs as ccrs #provides functionality for geographic projections
from cartopy.util import add_cyclic_point #provides functionality for closing gaps at periodic boundaries
from cdo import * #provides python bindings for climate data operators (CDO)
cdo=Cdo() #towards easier calling of Cdo()
#for further information regarding CDO see https://code.mpimet.mpg.de/projects/cdo/

defining various parameters to this script

#definition of filename templates

#for experiment data files
exp_file_name_template="data/MODEL_EPOCH_XxCO2_100years.nc" #place holders "MODEL", "EPOCH", and "X" will be replaced further down with simulation specific strings
#experiment files available:
#COSMOS-AO_Maastrichtian_2xCO2_100years.nc
#COSMOS-AO_Miocene_2xCO2_100years.nc
#COSMOS-AO_Oligocene_2xCO2_100years.nc
#COSMOS-AO_PreIndustrial_2xCO2_100years.nc
#COSMOS-ASO_Pliocene_2xCO2_100years.nc
#COSMOS-ASO_PreIndustrial_2xCO2_100years.nc
#MPI-ESM_PreIndustrial_2xCO2_100years.nc
#MPI-ESM_PreIndustrial_4xCO2_100years.nc

#for reference data files
ref_file_name_template="data/MODEL_PreIndustrial_100years.nc" #place holders "MODEL", "EPOCH", and "X" will be replaced further down with simulation specific strings
#reference files available:
#COSMOS-AO_PreIndustrial_100years.nc
#COSMOS-ASO_PreIndustrial_100years.nc
#MPI-ESM_PreIndustrial_100years.nc

#definition of parameters identifying experiment and reference data
model="COSMOS-AO"
epoch="Maastrichtian"
epoch_ref="PreIndustrial"
x="2"

setting up filenames based on above defined parameters

#generating filenames from filename templates and parameters
exp_file=exp_file_name_template.replace("MODEL",model).replace("EPOCH",epoch).replace("X",x)
ref_file=ref_file_name_template.replace("MODEL",model)

preparing data by computing various averages in time and space

We employ the Climate Data Operators (https://code.mpimet.mpg.de/projects/cdo/) that provide methods optimized for computing correctly area-weighted spatial averages of climate characteristics on a sphere. Furthermore, the CDO are capable of correctly computing time averages for the Gregorian Calendar, considering that different months have different lengths. The CDO are able to read, interpret, and write gridded binary climate data from the binary Network Common Data Format (netCDF, https://www.unidata.ucar.edu/software/netcdf/) that has been used to store the climate data at hand.

#computing various time and spatial averages for experiment and reference
filename_dictionary={
  "yearmeans": "",
  "timmean": "",
  "yearmeans_fldmean": "",
  "timmean_fldmean": "",
  "yearmeans_arctic_mean": "",
  "timmean_arctic_mean": "",
  "yearmeans_arctic_mean": "",
  "timmean_arctic_mean": "",
  "yearmeans_lowlat_mean": "",
  "timmean_lowlat_mean": "",
  "yearmeans_arctic_mean": "",
  "timmean_arctic_mean": "",
}
exp_files=filename_dictionary.copy()
ref_files=filename_dictionary.copy()
for realm in exp_file,ref_file:
    #computation of year means of spatially resolved monthly data
    outfile_yearmeans=realm.replace("data/","").replace('.nc','_yearmonmean.nc')
    retval=cdo.yearmonmean(input = realm, output = outfile_yearmeans)
    #computation of time means of spatially resolved yearly data
    outfile_timmean=outfile_yearmeans.replace("data/","").replace('.nc','_timmean.nc')
    retval=cdo.timmean(input = outfile_yearmeans, output = outfile_timmean)
    #computation of global mean of year average spatially resolved data
    outfile_yearmeans_fldmean=outfile_yearmeans.replace('.nc','_fldmean.nc')
    retval=cdo.fldmean(input = outfile_yearmeans, output = outfile_yearmeans_fldmean)
    #computation of global mean of time average spatially resolved data
    outfile_timmean_fldmean=outfile_timmean.replace('.nc','_fldmean.nc')
    retval=cdo.fldmean(input = outfile_timmean, output = outfile_timmean_fldmean)
    #computation of Arctic mean of year average spatially resolved data
    outfile_yearmeans_arctic_mean=outfile_yearmeans.replace('.nc','_arctic_mean.nc')
    retval=cdo.fldmean(input = "-sellonlatbox,0,360,60,90 %s"%(outfile_yearmeans), output = outfile_yearmeans_arctic_mean)
    #computation of Arctic mean of time average spatially resolved data
    outfile_timmean_arctic_mean=outfile_timmean.replace('.nc','_arctic_mean.nc')
    retval=cdo.fldmean(input = "-sellonlatbox,0,360,60,90 %s"%(outfile_timmean), output = outfile_timmean_arctic_mean)
    #computation of low latitude mean of year average spatially resolved data
    outfile_yearmeans_lowlat_mean=outfile_yearmeans.replace('.nc','_lowlat_mean.nc')
    retval=cdo.fldmean(input = "-sellonlatbox,0,360,0,30 %s"%(outfile_yearmeans), output = outfile_yearmeans_lowlat_mean)
    #computation of low latitude mean of time average spatially resolved data
    outfile_timmean_lowlat_mean=outfile_timmean.replace('.nc','_lowlat_mean.nc')
    retval=cdo.fldmean(input = "-sellonlatbox,0,360,0,30 %s"%(outfile_timmean), output = outfile_timmean_lowlat_mean)
    if realm == exp_file:
        exp_files["yearmeans"]=outfile_yearmeans
        exp_files["timmean"]=outfile_timmean
        exp_files["yearmeans_fldmean"]=outfile_yearmeans_fldmean
        exp_files["timmean_fldmean"]=outfile_timmean_fldmean
        exp_files["yearmeans_lowlat_mean"]=outfile_yearmeans_lowlat_mean
        exp_files["timmean_lowlat_mean"]=outfile_timmean_lowlat_mean
        exp_files["yearmeans_arctic_mean"]=outfile_yearmeans_arctic_mean
        exp_files["timmean_arctic_mean"]=outfile_timmean_arctic_mean
    if realm == ref_file:
        ref_files["yearmeans"]=outfile_yearmeans
        ref_files["timmean"]=outfile_timmean
        ref_files["yearmeans_fldmean"]=outfile_yearmeans_fldmean
        ref_files["timmean_fldmean"]=outfile_timmean_fldmean
        ref_files["yearmeans_lowlat_mean"]=outfile_yearmeans_lowlat_mean
        ref_files["timmean_lowlat_mean"]=outfile_timmean_lowlat_mean
        ref_files["yearmeans_arctic_mean"]=outfile_yearmeans_arctic_mean
        ref_files["timmean_arctic_mean"]=outfile_timmean_arctic_mean

loading data from disk

We employ the Dataset structure of the netCDF4 package that provides high level access to data stored in the NetCDF binary format.

#load time average spatially resolved data
exp_file_hdl=Dataset(exp_files["timmean"]) #open NetCDF file named in entry *timmean* of dictionary *exp_files*
ref_file_hdl=Dataset(ref_files["timmean"])
exp_temp2_timmean=exp_file_hdl.variables['temp2'][:].squeeze() #read specific variable from the opened NetCDF file
ref_temp2_timmean=ref_file_hdl.variables['temp2'][:].squeeze()
exp_file_hdl.close() #close the NetCDF file
ref_file_hdl.close()

#load yearmean global average time series
exp_file_hdl=Dataset(exp_files["yearmeans_fldmean"])
ref_file_hdl=Dataset(ref_files["yearmeans_fldmean"])
exp_temp2_yearmeans_fldmean=exp_file_hdl.variables['temp2'][:].squeeze()
ref_temp2_yearmeans_fldmean=ref_file_hdl.variables['temp2'][:].squeeze()
exp_file_hdl.close()
ref_file_hdl.close()

#load time average global average
exp_file_hdl=Dataset(exp_files["timmean_fldmean"])
ref_file_hdl=Dataset(ref_files["timmean_fldmean"])
exp_temp2_timmean_fldmean=exp_file_hdl.variables['temp2'][:].squeeze()
ref_temp2_timmean_fldmean=ref_file_hdl.variables['temp2'][:].squeeze()
exp_file_hdl.close()
ref_file_hdl.close()

#load yearmean low latitude average
exp_file_hdl=Dataset(exp_files["yearmeans_lowlat_mean"])
ref_file_hdl=Dataset(ref_files["yearmeans_lowlat_mean"])
exp_temp2_yearmeans_lowlat_mean=exp_file_hdl.variables['temp2'][:].squeeze()
ref_temp2_yearmeans_lowlat_mean=ref_file_hdl.variables['temp2'][:].squeeze()
exp_file_hdl.close()
ref_file_hdl.close()

#load time average low latitude average
exp_file_hdl=Dataset(exp_files["timmean_lowlat_mean"])
ref_file_hdl=Dataset(ref_files["timmean_lowlat_mean"])
exp_temp2_timmean_lowlat_mean=exp_file_hdl.variables['temp2'][:].squeeze()
ref_temp2_timmean_lowlat_mean=ref_file_hdl.variables['temp2'][:].squeeze()
exp_file_hdl.close()
ref_file_hdl.close()

#load yearmean arctic average
exp_file_hdl=Dataset(exp_files["yearmeans_arctic_mean"])
ref_file_hdl=Dataset(ref_files["yearmeans_arctic_mean"])
exp_temp2_yearmeans_arctic_mean=exp_file_hdl.variables['temp2'][:].squeeze()
ref_temp2_yearmeans_arctic_mean=ref_file_hdl.variables['temp2'][:].squeeze()
exp_file_hdl.close()
ref_file_hdl.close()

#load time average arctic average
exp_file_hdl=Dataset(exp_files["timmean_arctic_mean"])
ref_file_hdl=Dataset(ref_files["timmean_arctic_mean"])
exp_temp2_timmean_arctic_mean=exp_file_hdl.variables['temp2'][:].squeeze()
ref_temp2_timmean_arctic_mean=ref_file_hdl.variables['temp2'][:].squeeze()
exp_file_hdl.close()
ref_file_hdl.close()

#load geographic coordinates (for plotting geographic maps)
#it is assumed that reference data and experiment data are on the same grid, hence coordinates are loaded only once
exp_file_hdl=Dataset(exp_files["timmean"])
ref_file_hdl=Dataset(ref_files["timmean"])
model_lon=exp_file_hdl.variables['lon'][:]
model_lat=exp_file_hdl.variables['lat'][:]
exp_file_hdl.close()
ref_file_hdl.close()

plotting the climate state using various aggregation methods

global map of multiyear time average anomaly of near surface air temperature

#plot geograpic map of anomaly exp vs. ref
temp2_anomaly=exp_temp2_timmean-ref_temp2_timmean
temp2_anomaly_cyclic, model_lon_cyclic=add_cyclic_point(temp2_anomaly, coord=model_lon) #close the data gap by adding cyclic point
fig, ax = plt.subplots(1, 1, dpi=300, subplot_kw={"projection": ccrs.Robinson(central_longitude=0)})
ax.coastlines(linewidth=0.33,color='gray')
ax.set_global()
ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
             linewidth=0.5, color='gray', alpha=0.5, linestyle='-')
level_min=-15
level_max=15
stepsize=1
levels=np.arange(level_min-stepsize/2,level_max+stepsize/2+stepsize,stepsize)
cmap=colors.LinearSegmentedColormap.from_list('mycmap', 
             ['blue', 'cornflowerblue', 'aqua', 'lightcyan', 'white', 'yellow', 'orange', 'red', 'maroon'],
             (level_max-level_min)/stepsize+1) #+1 adds one level for the center
cf=ax.contourf(model_lon_cyclic, model_lat, temp2_anomaly_cyclic, transform=ccrs.PlateCarree(),
               cmap=cmap, levels=levels, 
               vmin=level_min-.5, vmax=level_max+.5,
               extend='both')
cax=plt.colorbar(cf, ticks=np.arange(level_min+1,level_max+1-1,2), drawedges=True, orientation='horizontal',
                       label=r"near surface air temperature anom. / K",pad=0.05,shrink=1.25,aspect=40)
cax.ax.xaxis.set_ticklabels([label.get_text().replace("−", "-") for label in cax.ax.xaxis.get_ticklabels()])
cax.ax.tick_params(labelsize=9)

Fig. 4: Anomaly of near surface air temperature, experiment minus reference, based on a multiyear annual mean computed over 100 model years.

annual mean time series of absolute values of near air surface temperature, spatially averaged in the global domain

#plot yearmean global mean time series for experiment and reference
plt.plot(exp_temp2_yearmeans_fldmean,marker='o',linewidth=3,label=epoch,color="C1")
plt.plot(np.repeat(np.mean(exp_temp2_yearmeans_fldmean),len(exp_temp2_yearmeans_fldmean)),linewidth=3,label="time average "+epoch+" ("+str(np.round(np.mean(exp_temp2_yearmeans_fldmean),2))+" K)",alpha=0.8,color="C1")
plt.plot(ref_temp2_yearmeans_fldmean,marker='o',linewidth=3,label=epoch_ref,color="C0")
plt.plot(np.repeat(np.mean(ref_temp2_yearmeans_fldmean),len(ref_temp2_yearmeans_fldmean)),linewidth=3,label="time average "+epoch_ref+" ("+str(np.round(np.mean(ref_temp2_yearmeans_fldmean),2))+" K)",alpha=0.8,color="C0")
plt.xlabel("year")
plt.ylabel("near surface air temperature / K")
plt.title("yearmean global average")
plt.legend()
plt.show()

Fig. 5: Absolute values of annual average global average near surface air temperature for both experiment and reference. Shown is the available time series that covers 100 model years.

annual mean time series of absolute values of near surface air temperature, spatially averaged across the Arctic

#plot yearmean Arctic mean time series for experiment and reference
plt.plot(exp_temp2_yearmeans_arctic_mean,marker='o',linewidth=3,label=epoch,color="C1")
plt.plot(np.repeat(np.mean(exp_temp2_yearmeans_arctic_mean),len(exp_temp2_yearmeans_arctic_mean)),linewidth=3,label="time average "+epoch+" ("+str(np.round(np.mean(exp_temp2_yearmeans_arctic_mean),2))+" K)",alpha=0.8,color="C1")
plt.plot(ref_temp2_yearmeans_arctic_mean,marker='o',linewidth=3,label=epoch_ref,color="C0")
plt.plot(np.repeat(np.mean(ref_temp2_yearmeans_arctic_mean),len(ref_temp2_yearmeans_arctic_mean)),linewidth=3,label="time average "+epoch_ref+" ("+str(np.round(np.mean(ref_temp2_yearmeans_arctic_mean),2))+" K)",alpha=0.8,color="C0")
plt.xlabel("year")
plt.ylabel("near surface air temperature / K")
plt.title("yearmean Arctic average")
plt.legend()
plt.show()

Fig. 6: As Fig. 5 but for the Arctic average.

annual mean time series of absolute values of near surface air temperature, spatially averaged across the low latitudes of the Northern Hemisphere

#plot yearmean low latitude mean time series for experiment and reference
plt.plot(exp_temp2_yearmeans_lowlat_mean,marker='o',linewidth=3,label=epoch,color="C1")
plt.plot(np.repeat(np.mean(exp_temp2_yearmeans_lowlat_mean),len(exp_temp2_yearmeans_lowlat_mean)),linewidth=3,label="time average "+epoch+" ("+str(np.round(np.mean(exp_temp2_yearmeans_lowlat_mean),2))+" K)",alpha=0.8,color="C1")
plt.plot(ref_temp2_yearmeans_lowlat_mean,marker='o',linewidth=3,label=epoch_ref,color="C0")
plt.plot(np.repeat(np.mean(ref_temp2_yearmeans_lowlat_mean),len(exp_temp2_yearmeans_lowlat_mean)),linewidth=3,label="time average "+epoch_ref+" ("+str(np.round(np.mean(ref_temp2_yearmeans_lowlat_mean),2))+" K)",alpha=0.8,color="C0")
plt.xlabel("year")
plt.ylabel("near surface air temperature / K")
plt.title("yearmean low latitude average")
plt.legend()
plt.show()

Fig. 7: As Fig. 5 but for the Northern Hemisphere low-latitude average.

derive climate metrics from quantities computed above

We compute here two climate metrics that refer to important characteristics of the Earth's climate:

• Climate Sensitivity (or Earth System Sensitivity in case experiment and reference are based on different sets of land surface conditions), that is defined as the global average temperature change occurring due to a doubling of carbon dioxide from 280 ppmv to 560 ppmv
• Arctic Amplification, that is defined as the factor by which the Arctic warms more than the low latitude regions of the Earth do

#climate sensitivity (global average temperature change for a doubling of carbon dioxide)
print("Climate / Earth System sensitivity for epoch "+epoch+", with reference to epoch "+epoch_ref+": "+str(np.round(exp_temp2_timmean_fldmean-ref_temp2_timmean_fldmean,2))+" K")

Climate / Earth System sensitivity for epoch Maastrichtian, with reference to epoch PreIndustrial: 5.7 K

#Arctic amplification (difference in warming at northern high latitudes vs. northern low latitudes)
print("Arctic Amplification for epoch "+epoch+", with reference to epoch "+epoch_ref+": "+str(np.round((exp_temp2_timmean_arctic_mean-ref_temp2_timmean_arctic_mean)/(exp_temp2_timmean_lowlat_mean-ref_temp2_timmean_lowlat_mean),2))+" K")

Arctic Amplification for epoch Maastrichtian, with reference to epoch PreIndustrial: 2.42 K

Further tasks

You may extend analysis of the climate states by performing some or all of the following tasks:

• Analyse the impact of a doubling carbon dioxide on global precipitation (variable "precip" in files data/{COSMOS-AO_PreIndustrial_100years.nc,COSMOS-AO_Maastrichtian_2xCO2_100years.nc}).
• In case you analyzed precipitation as suggested above: By which fraction (as a percentage of the precipitation in the reference state) does precipitation increase per degree centigrade of global average warming?
• Estimate uncertainty in the values derived for Climate Sensitivity (Earth System Sensitivity) and Polar Amplification; employ any of the following methods:
  • propagate uncertainties in "measured" quantities (e.g. standard deviation over time of global average temperatures in experiment and reference) to derive uncertainty in computed quantities (Gaussian Error Propagation)
  • derive uncertainty by computing the relevant climate metric from annual values (rather than from the multiyear annual mean over 100 years as performed in this notebook) and derive standard deviation from the resulting time series of Climate Sensitivity, Earth System Sensitivity, or Arctic Amplification, respectively.
• Quantify the relationship between Climate Sensitivity, respectively Earth System Sensitivity, and Polar Amplification. To this end plot tuples of both quantities for various simulations in the phase space and derive their relationship via regression analysis - ask your peers for values of Climate Sensitivity / Earth System Sensitivity and Polar Amplification they derived from their simulations)
• How large is Polar Amplification in the Southern Hemisphere?