Unrealistic oxygen isotopic composition of precipitation in iCESM1.2?

[AndreaMoore]

Hello,

I recently completed some paleoclimate simulations in iCESM1.2, and am now trying to compute d18O_precipitation values from the output using the following code:

# Convert isotope data from rain rate (in m/s) to d18O
fill_val = np.nan   # old fill_val = 9.96921e+36
Rstd = 2005.20e-6   # O18/O16 ratio in standard (here, VSMOW) = 2005.20 ± 0.45 ppm (Sharp txtbook, 2007)
Mh2o = 18.015       # mean molecular mass of water
Mh218o = 20.013     # molecular mass of H218O, units = atomic mass units (1/12 mass of 12C ~ 1.66654 E -27 kg)
Mh216o = 18.009     # molecular mass of H216O

# Load in variables as values
H216ORL_ctrl = ds_H216ORL_ctrl.PRECRL_H216OR.values
H216ORC_ctrl = ds_H216ORC_ctrl.PRECRC_H216Or.values
H218ORL_ctrl = ds_H218ORL_ctrl.PRECRL_H218OR.values
H218ORC_ctrl = ds_H218ORC_ctrl.PRECRC_H218Or.values

# Add large-scale and convective values to get total
H216O_ctrl = H216ORL_ctrl + H216ORC_ctrl
H218O_ctrl = H218ORL_ctrl + H218ORC_ctrl
H216O_ctrl[H216O_ctrl == 0] = fill_val # Replace 0's with NaNs for R18O calculation

R18O_ctrl = (H218O_ctrl/H216O_ctrl)*(Mh216o/Mh218o)
d18O_ctrl = (R18O_ctrl - 1)*1000 # Rstandard assumed = 1

However, I'm getting some unrealistic d18O_precip values out of this calculation, ranging from -9148167.0 permil to 77604790.0 permil. Am I doing something wrong? Should I simply set a threshold value and ignore anything above and below it?

Any help is appreciated!

 

[nusbaume]

Hi Andrea,

There are few things I would do to see if I can get reasonable-looking isotope values.

First, I would make sure that the raw “H216O” and “H218O” fields look reasonable, i.e. convert them to mm/day and just see if a plot of them looks like what you would expect regular model rain to also look like (you can compare them to "PRECC" and "PRECL" if you want, although note those also include “snow” values as well, which can be removed using "PRECSC" and "PRECSL" if need be).

If they look reasonable, then the next thing I would do is remove all isotope values before summing if their mass total is below a small, but non-zero, threshold, i.e. something like this:

H216ORL_ctrl[H216ORL_ctrl < 10**-18] = 0
H216ORC_ctrl[H216ORC_ctrl < 10**-18] = 0
H218ORL_ctrl[H218ORL_ctrl < 10**-18] = 0
H2168ORC_ctrl[H218ORC_ctrl < 10**-18] = 0


H216O_ctrl = H216ORL_ctrl + H216ORC_ctrl
H218O_ctrl = H218ORL_ctrl + H218ORC_ctrl

Then you’ll want to do the same thing for the sum like you are already doing, except you’ll also want a small non-zero number as the threshold:

H216O_ctrl[H216O_ctrl < 10**-18] = fill_val

The reason for all of this is that the model numerics can produce unrealistic isotope values at very small (but still non-zero) mass values, which is what all of this post-processing code is trying to avoid.

Finally, you don’t want to multiply with the molecular weight ratio, as that is also “included” in the "Rstandard = 1" conversion.

Anyways, I hope that helps, and if you still are getting bad values after these adjustments then please let me know.

Thanks, and have a great day!

Jesse

 

[AndreaMoore]

Hi Jesse,

Thanks for the help! Everything seems to be working now.

Cheers,
Andrea