New recipe: rate of sea ice area loss per degree warming #3983
Conversation
Only partial functionality at this stage
|
Dear @alistairsellar , thank you for starting on this. I am the PI for the OSISAF/CCI sea-ice concentration/area data and have previously assisted @axel-lauer for the other sea-ice related metric (seasonal cycles of SIA). I am happy to also try and assist in this metric. Is it the plan to use the same timeserie(s) of observed sea-ice area as for the other metric? |
Hi @TomLav thanks for your interest and offer of help! Would you be willing to act as a science reviewer?
For this PR we are planning to hard-code uncertainty ranges from the papers we are reproducing - the ranges are based on a combination of obs and model data. See #3962 (comment). If you were keen to get the observations into a subsequent implementation of the diagnostic, you would be welcome to explore reproducing the uncertainty calculation of the papers listed in #3962. The result would give a different answer from the papers since the data would be different, but the recipe could include an option to switch between the hard-coded values from the original papers and an updated value calculated from the latest obs data and models selected by the user. |
I think this is a fair approach. Do you have access to the underlying data the authors used for the two CMIP6 papers? They shared their software (data?) for the Arctic paper, at least: https://github.com/jakobdoerr/SIMIP_2020/tree/master
An update might come later, indeed. From what I see, the metric in the papers cover up-until 2014 (historical period?), and I do not know if this is still ok for CMIP7.
My expertise is really only on the sea-ice observation side of things, so I am not sure I am a good reviewer for this sensitivity metric. I am also not a user of ESMValTool. Maybe you can hear if other (sea ice) climate modellers would be available as reviewers and I would take the job if none have the time? |
This is really part of the last commit but my SSH to GitHub seems to be broken
logging is a guess, provenance record is not accurate
As this seems to be used in the Roach et al plot
I have queries about this to send to Alistair
This wasn't meant to be meaned globally
There was a problem hiding this comment.
Hi @NParsonsMO, the code is looking great! A couple of suggestions and notes in these comments.
| logger.debug(f'Calculating linear relationship between {dependent} and {independent}') | ||
|
|
||
| # Use SciPy stats to calculate the regression | ||
| slope, intercept, r, p, stderr = stats.linregress(independent, dependent) |
There was a problem hiding this comment.
We might consider using the p-value from this to reproduce the hatching capability of Roach's Figure 3e.
There was a problem hiding this comment.
Added in latest version.
It seemed silly to have three different functions all just calling stats.linregress so I replaced them with just one.
| x = 10 * inner_dict['tas_trend'] # x10 to approximate decades | ||
| y = 10 * inner_dict['siconc_trend'] # x10 to approximate decades |
There was a problem hiding this comment.
Can we use cube.convert_units("K decade-1") or a convert_units preprocessor instead of manual conversion?
There was a problem hiding this comment.
I am struggling to understand, sorry...
The cube's units were just K (or area, in a couple of different forms, for siconc). By the time there was a trend (K/year) then we only had a single numerical value for slope (calculated with linregress) rather than a cube.
Do you want the cube to be transformed into having a smaller length in the time dimension, but with (presumably averaged) data points to represent an entire decade? [I was under the impression, "no."]
|
Hi @ESMValGroup/tech-reviewers please could someone volunteer or suggest a technical reviewer for this PR, which provides a diagnostic required for the REF? It's still in draft but will be ready for review soon. |
I can have a look at it later! In the meantime, could you please run |
|
Hi @schlunma, Thanks for the offer to review! I have a couple of changes left to make to the plotting functions* but otherwise I think it is hopefully ready to review, and they shouldn't really affect anything... I won't be able to work on this for a couple of weeks, so it might be worth looking beforehand! Not least because there might well be other stuff that I will then fix at the same time. Thanks again!
|
There was a problem hiding this comment.
Thanks @NParsonsMO and @alistairsellar for your contribution! Code looks great, just have a couple of minor comments.
This is ready for scientific review now! 🚀
Note: I did not run the recipe.
There was a problem hiding this comment.
Would it be possible to ensure that the title is not cut off? This can probably be achieved by passing bbox_inches="tight" as argument for plt.savefig.
It would also be helpful to provide legend entries for the dashed lines and shaded region.
|
|
||
| * mean: -4.01, | ||
| * std_dev: 0.32, | ||
| * plausible: 1.28, |
There was a problem hiding this comment.
It's probably better to avoid bullet points here, these are rendered very awkwardly.
| assert len(selection) == 1, ( | ||
| f"None or too many matching files found for {variable_group}" | ||
| ) |
There was a problem hiding this comment.
Instead of assert, please use exceptions (e.g., a raise ValueError) to validate user input.
| defaults: &defaults {ensemble: r1i1p1f1, exp: historical, grid: gn, project: CMIP6} | ||
| datasets: | ||
| - {<<: *defaults, dataset: HadGEM3-GC31-LL, institute: MOHC, ensemble: r1i1p1f3} | ||
| # - {<<: *defaults, dataset: HadGEM3-GC31-MM, institute: MOHC, ensemble: r1i1p1f3} # other error |
There was a problem hiding this comment.
Please remove all commented lines.
schlunma
added
looking for scientific reviewer
approved by technical reviewer
diagnostic
new recipe
| "notz_fig_title": "September Arctic Sea Ice Sensitivity", | ||
| "notz_ax_title": f"dSIA/dGMST obs from {next(iter(dictionary['obs'].keys()))}", |
There was a problem hiding this comment.
| "notz_fig_title": "September Arctic Sea Ice Sensitivity", | |
| "notz_ax_title": f"dSIA/dGMST obs from {next(iter(dictionary['obs'].keys()))}", | |
| "notz_fig_title": "September Arctic sea-ice area sensitivity to global mean surface temperature", | |
| "notz_ax_title": f"Annual regression over {next(iter(dictionary['obs'].keys()))}", |
Suggested more verbose titles.
I removed "obs" - is that OK or does it obscure an intended meaning?
There was a problem hiding this comment.
I think yours is better, but only the (currently hardcoded) observations are definitely from those years, the years of the data are configurable in the recipe (although they are currently set to the same), so I probably need a different title.
There was a problem hiding this comment.
After discussion: I plan to remove this axis title from here and put a more descriptive caption in place, including both the data period and observations period.
| region = [item["diagnostic"] for item in cfg["input_data"].values()][0] | ||
|
|
||
| record = { | ||
| "caption": f"Plots of {region.title()} sea ice sensitivity", |
There was a problem hiding this comment.
From discussion at ESMValTool workshop today there is a request for all REF-focused recipes to ensure that provenance captions are meaningful and specific, because the captions will be used in web publication of plots.
We therefore probably should have a different caption for each plot. I suggest that we add a caption argument to this function to allow the calling routines to specify it themselves. I recommend using the fig_title variables (without the ax_title subheadings.
There was a problem hiding this comment.
This is a science review ahead of finalising the code.
The code looks neat and the desired figures seem to be replicated reasonably well.
I didn't see any major issues but have left some comments on bits to change and asked for clarification on some possible discrepancies values in the Southern Ocean figure.
| :align: center | ||
| :width: 18cm | ||
|
|
||
| Plot of the trend of annually averaged southern hemisphere sea ice area (millions of square kilometres) over time against the trend of annually and globally averaged air temperature near the surface (degrees Kelvin) over time. |
There was a problem hiding this comment.
This line should include time ranges for the trends. Also, how ensemble members are handled should be described as that looks like a difference between this figure and Roach et al.
|
|
||
| Plot of the trend of annually averaged southern hemisphere sea ice area (millions of square kilometres) over time against the trend of annually and globally averaged air temperature near the surface (degrees Kelvin) over time. | ||
|
|
||
| The colour of each point is determined by the ``r value`` of the correlation between the two variables, and the hatching indicates a ``p value`` greater than ``0.05``, both calculated using ``scipy.stats.linregress``. |
There was a problem hiding this comment.
Here, r should be decribed as the Pearson correlation coefficient for precision and consistency with Roach et al.
There was a problem hiding this comment.
I think something is off with the title here. Also, the cbar caption should be reformatted to have an italicised R, superscript 2, and a multiplication sign rather than *. You may also want to consider placing the model labels to the left as there's more space there and they wouldn't push the cbar out.
The figure looks generally consistent with Roach et al, though there appear to be slightly more positive SIA trends in this figure. Is there a know reason for that? (see, for example, CAMS-CSM1-0 up around 0.2, where nothing is up to 0.1 in Roach et al.)
Is there a reason no obs are included?
There was a problem hiding this comment.
I think 'obs' should be removed from the title. Also either 'sea ice sensitivity' or 'dSIA/dGMST' could potentially be removed. Probably 'dSIA/dGMST' should be removed as it appears in the ylabel too.
Also, the observation mean and shading should be added to the legend. The dotted grey lines seem a bit faint (I didn't notice them on first look).
The figure looks generally consistent with Notz et al.
| :align: center | ||
| :width: 8cm | ||
|
|
||
| Plot of northern hemisphere sea ice area loss (millions of square kilometres) in the month of September per degree Kelvin of average global warming. |
There was a problem hiding this comment.
I would suggest replacing [per degree Kelvin of average global warming] with [per global mean temperature change (K)].
In general this text is a bit technical, and introducing the caption with something like 'Sensitivity of Arctic sea ice loss to...' would be helpful.
Also please add some mention of how ensemble members are handled.
|
|
||
| The dashed black line shows the observational mean, the shaded area is within one standard deviation of the mean, and the dashed grey lines are within "plausibility" of the mean, with all values taken from Ed Blockley's code for the period 1979-2014 and defined as follows: | ||
|
|
||
| * mean: -4.01, |
There was a problem hiding this comment.
These values are missing units.
Also added more logger info
Only partial functionality at this stage
Description
Implements a diagnostic and recipe for quantifying the sensitivity of seaice to global warming.
Before you get started
Checklist
It is the responsibility of the author to make sure the pull request is ready to review. The icons indicate whether the item will be subject to the 🛠 Technical or 🧪 Scientific review.
New or updated recipe/diagnostic
To help with the number of pull requests: