3: Water isotope tracers in CESM

3: Water isotope tracers in CESM#

Exercise: Run a preindustrial simulation with water isotope tracers

Download (git clone) isotope-enabled CESM1.3 (iCESM1.3; code available here), as the version of CESM used in this tutorial does not include water isotope capabilities.

Create, configure, build and run a fully coupled preindustrial case called b.e13.B1850.f19_g16.piControl.001 following CESM naming conventions including water isotope tracers.

Run for 1 year.

Click here for hints

What is the resolution for B1850?

Use resolution f19_g16 for fast throughput

Which XML variable should you change to tell the model to run for one year?

Use STOP_OPTION and STOP_N

How to check if each XML variable is modified correctly?

Use xmlquery -p

Click here for the solution

# Download iCESM1.3 code

cd /glade/u/home/$USER/code 
git clone https://github.com/NCAR/iCESM1.3_iHESP_hires iCESM1.3_iHESP_hires 
cd iCESM1.3_iHESP_hires 
./manage_externals/checkout_externals 

# Set environment variables

Set environment variables with the commands:

For tcsh users

set CASENAME=b.e13.B1850C5.f19_g16.piControl.001
set CASEDIR=/glade/u/home/$USER/cases/$CASENAME
set RUNDIR=/glade/derecho/scratch/$USER/$CASENAME/run
set COMPSET=B1850C5
set RESOLUTION=f19_g16
set PROJECT=UESM0013

Note: You should use the project number given for this tutorial.

For bash users

export CASENAME=b.e13.B1850C5.f19_g16.piControl.001
export CASEDIR=/glade/u/home/$USER/cases/$CASENAME
export RUNDIR=/glade/derecho/scratch/$USER/$CASENAME/run
export COMPSET=B1850C5
export RESOLUTION=f19_g16
export PROJECT=UESM0013

Note: You should use the project number given for this tutorial.

# Make a case directory

If needed create a directory cases into your home directory:

mkdir /glade/u/home/$USER/cases/

# Create a new case

Create a new case with the command create_newcase:

cd /glade/u/home/$USER/code/iCESM1.3_iHESP_hires/cime/scripts/
./create_newcase --case $CASEDIR --res $RESOLUTION --compset $COMPSET --project $PROJECT --run-unsupported 

# Change the job queue

If needed, change job queue.
For instance, to run in the queue main.

cd $CASEDIR
./xmlchange JOB_QUEUE=main

This step can be redone at anytime in the process before running case.submit.

# Setup

Invoke case.setup with the command:

cd $CASEDIR
./case.setup    

You build the namelists with the command:

./preview_namelists

This step is optional as the script preview_namelists is automatically called by case.build and case.submit. But it is nice to check that your changes made their way into:

$CASEDIR/CaseDocs/atm_in

# Set run length

./xmlchange STOP_N=1,STOP_OPTION=nyears

# Build the run

qcmd -A $PROJECT -- ./case.build

# Which namelist variables enable water isotope tracers?

Notice that the steps to set up this isotope-enabled preindustrial simulation are very similar to a preindustrial simulation without isotopes (e.g., Paleo Exercise 1)
In iCESM1.3, it is assumed you will run with water isotope tracers so each compset include isotope settings by default
Use ./xmlquery to explore how FLDS_WISO, CAM_CONFIG_OPTS, and OCN_TRACER_MODULES differ between this iCESM1.3 case and that of Exercise 1
Also, take a look at namelist settings for each CESM component with variables that contain wiso in $CASEDIR/Buildconf

# Submit the run

./case.submit

# Check on your run

After submitting the job, use qstat -u $USER to check the status of your job.

# Check your solution

When the run is completed, look at the history files into the archive directory.

(1) Check that your archive directory on derecho (The path will be different on other machines):

cd /glade/derecho/scratch/$USER/archive/$CASENAME/atm/hist
ls 

As your run is one-year, there should be 12 monthly files (h0) for each model component.

Success! Let’s plot the results.

Click here to visualize results

————— Option 1 —————

# Use NCO to calculate the oxygen isotopic composition of precipitation

The ratio of heavy ($^{18}\text{O}$) to light ($^{16}\text{O}$)) isotopes are most commonly expressed relative to a standard in delta (δ) notation:

\[ \delta^{18}\text{O} = \frac{R_{\text{sample}} - R_{\text{std}}}{R_{\text{std}}} \times 1000‰ \]

where

$R_{\text{sample}}$ = ratio of $^{18}\text{O}$ to $^{16}\text{O}$ in sample
$R_{\text{std}}$ = ratio of $^{18}\text{O}$ to $^{16}\text{O}$ in a standard

Thus, the $\delta^{18}\text{O}$ of a sample which is identical to the standard would be 0‰, positive values indicate a greater proportion of $^{18}\text{O}$ than the standard, and negative values indicate a lower proportion of $^{18}\text{O}$ .

In isotope-enabled CESM, the relative abundances of $^{16}\text{O}$ and $^{18}\text{O}$ are already adjusted to their naturally occurring global abundances (99.757% and 0.205%, respectively), so we do not include $R_{\text{std}}$ in the calculation of $\delta^{18}\text{O}$. Rather, isotope variables in CESM are expressed in delta (δ) notation as:

\[ \delta^{18}O = (\frac{\text{PRECRC_H218Or} + \text{PRECSC_H218Os} + \text{PRECRL_H218OR} + \text{PRECSL_H218OS}}{\text{PRECRC_H216Or} + \text{PRECSC_H216Os} + \text{PRECRL_H216OR} + \text{PRECSL_H216OS}} - 1) \times 1000‰ \]

Use ncdump /glade/derecho/scratch/$USER/$CASENAME/atm/hist/$CASENAME.cam.h0.0001-01.nc | less to check the definition of each isotope variable above
For example, search for PRECRC_H218Or using /PRECRC_H218Or + <enter>)

To calculate the δ¹⁸O of precipitation from the simulation using NCO,

cd /glade/derecho/scratch/$USER/archive/$CASENAME/atm/hist
ncap2 -s 'd18Op=((PRECRC_H218Or+PRECSC_H218Os+PRECRL_H218OR+PRECSL_H218OS)/(PRECRC_H216Or+PRECSC_H216Os+PRECRL_H216OR+PRECSL_H216OS) - 1)*1000.' -v $CASENAME.cam.h0.0001-12.nc d18Op.$CASENAME.cam.h0.0001-12.nc 

# Use Ncview to visualize precipitation $\delta^{18}\text{O}$

Earth’s orbital configuration influences incoming solar insolation. Take a look at the d18Op variable we calculated for 1 month in the pre-industrial run.

module load ncview
ncview d18Op.$CASENAME.cam.h0.0001-12.nc 

————— Option 2 —————

# Use Python to calculate and plot the oxygen isotopic composition of precipitation

The following Python code will produce a plot of precipitation δ¹⁸O for 1 month.

import xarray as xr 
import numpy as np 
import matplotlib.pyplot as plt 
import cartopy.crs as ccrs 
from cartopy.util import add_cyclic_point 

def calculate_d18Op(ds): 
    # Compute precipitation δ18O with iCESM output 
    
    # Parameters 
    # ds: xarray.Dataset contains necessary variables 
    
    # Returns 
    # ds: xarray.Dataset with δ18O added 
    
    # convective & large-scale rain and snow, respectively 
    p16O = ds.PRECRC_H216Or + ds.PRECSC_H216Os + ds.PRECRL_H216OR + ds.PRECSL_H216OS 
    p18O = ds.PRECRC_H218Or + ds.PRECSC_H218Os + ds.PRECRL_H218OR + ds.PRECSL_H218OS 
    
    # avoid dividing by small number here 
    p18O = p18O.where(p16O > 1.E-18, 1.E-18) 
    p16O = p16O.where(p16O > 1.E-18, 1.E-18) 
    d18O = (p18O / p16O - 1.0) * 1000.0 
    
    ds['p16O'] = p16O 
    ds['p18O'] = p18O 
    ds['d18O'] = d18O 
    return ds 

# Read in monthly file 
case = 'b.e13.B1850C5.f19_g16.piControl.001' 
file = '/glade/derecho/scratch/macarew/archive/'+case+'/atm/hist/'+case+'.cam.h0.0001-12.nc' 
ds = xr.open_mfdataset(file, parallel=True, 
                       data_vars='minimal', 
                       coords='minimal', 
                       compat='override') 

# Call function to add preciptation d18O to dataset 
ds = calculate_d18Op(ds) 

fig, ax = plt.subplots( 
    nrows=1, ncols=1, 
    figsize=(6, 2), 
    subplot_kw={'projection': ccrs.Robinson(central_longitude=210)}, 
    constrained_layout=True) 

d18O_new, lon_new = add_cyclic_point(ds.d18O[0,:,:], ds.lon) 

# Plot model results using contourf 
p0 = ax.contourf(lon_new, ds.lat, d18O_new, 
                 levels=np.linspace(-30, 0, 16), extend='both', 
                 transform=ccrs.PlateCarree()) 

plt.colorbar(p0, ax=ax) 
ax.set_title('Dec δ18Op of PI') 
ax.coastlines(linewidth=0.5) 

————— Questions for reflection —————

Do you notice any spatial patterns in precipitation $\delta^{18}\text{O}$?