CMIP7 historical stratospheric aerosol optical properties and volcanic stratospheric sulfur emissions

[thomasaubry]

This discussion is the place for all things stratospheric aerosol optical properties and volcanic stratospheric sulfur emissions! Version 2.2.1 of our dataset is on ESGF and is the version recommended for use in CMIP7 Fast Track. After an update in January 2026, CMIP7 ScenarioMIP runs should use version 2.2.2 of our Scenario files.

The aerosol optical property dataset is documented in (please contribute to the open discussion in GMD):
Aubry, T. J., Toohey, M., Khanal, S., Chim, M. M., Verkerk, M., Johnson, B., Schmidt, A., Kovilakam, M., Sigl, M., Nicholls, Z., Thomason, L., Naik, V., Rieger, L., Stiller, D., Ziegler, E., and Smith, I.: Stratospheric aerosol forcing for CMIP7 (part 1): Optical properties for pre-industrial, historical, and scenario simulations (version 2.2.1), EGUsphere [preprint], https://doi.org/10.5194/egusphere-2025-4990, 2025.

The volcanic emission dataset is documented in (please contribute to the open discussion in GMD):
Aubry, T. J., Sigl, M., Toohey, M., Chim, M. M., Verkerk, M., Schmidt, A., and Carn, S. A.: Stratospheric aerosol forcing for CMIP7 (part 2): Volcanic sulfur dioxide emissions, EGUsphere [preprint], https://doi.org/10.5194/egusphere-2026-546, 2026.

All codes and source datasets used to produce our datasets are available on zenodo and GitHub (minus GloSSAC which was too big to upload on GitHub). We are keeping our brief documentation online where we will continue to document new versions of the dataset. Should any comment/question you leave be left unanswered for longer than you wish, please contact Thomas Aubry (thomas.aubry@earth.ox.ac.uk).

[MartineMichou]

Dear Dr. Aubry,

Dealing with the stratospheric volcano information you provided for CMIP6Plus we have two questions:

• in the CMIP6 volcano dataset there was an H2SO4 _mass field (molecules/cm3air) that corresponded to H2SO4 in the condensed phase.
  As our climate model uses this variable can you provide it on the same coordinates as the other fields?

• again, in the CMIP6 volcano dataset there were two sad fields, one called sad, and the other one sad_of_big_particles. We used the sad_of_big_particles field in our CMIP6 simulations.
  Can you provide sad_of_big_particles in addition to sad?

Looking forward to your response,

Martine Michou

[thomasaubry]

Dear Martine,
Replying on Github as well so that anyone with similar questions can see the answers.

1. For the aerosol number density (former H2SO4_msss field), we will provide it in the next version. I'm hoping to share an update through my own archive by the end of this week (Feb 7) and it should be on ESGF by early next week (~Feb 11). If you need me to accelerate this timeline please let me know and I will do my best (but travelling/working part time this week)
2. For the sad_of_big_particles field, it sounded from our discussion that you might not need it in the end. SHould any modelling group need this variable, please point me to documentation that enables me to understand your need as I cannot find this variable in any CMIP6 documentation.
   Merci!
   Thomas

[thomasaubry]

Hello,

Just confirming the number density is now added in version 1.3.1. This version is undergoing validation before being uploaded to ESGF, it should be there by Weds Feb 12 I suspect. If anyone wants a copy of the files (which at the moment have not been quality checked), please email t.aubry@exeter.ac.uk and I should come back to you on the same day.

All the best,

Thomas

[mdeushiMRI]

Dear Dr. Aubry,

Our climate model also uses the H2SO4 _mass field (molecules/cm3air).
Could you provide H2SO4 _mass data, please?

Thank you in advance.

Makoto Deushi.

[thomasaubry]

Dear Makoto,
Thanks a lot for this note. Please see my reply to Martine and let me know if you have any concern with the timeline. If you provide me with your email address (email t.aubry@exeter.ac.uk) I can share updated files before they are quality-checked and uploaded to ESGF.
All the best,
Thomas

[mdeushiMRI]

Dear Dr. Aubry,

Thank you for letting us know.

Just confirming the number density is now added in version 1.3.1.

We need (not the number density of ambient aerosol particles but) the H2SO4_mass field.
In the CMIP6 aerosolProperties data (ex. multiple_input4MIPs_aerosolProperties_CMIP_IACETH-SAGE3lambda-3-0-0_gn_185001_201412.nc) , the definition of the H2SO4_mass field is as follows.

       float H2SO4_mass(latitude, altitude, time) ;
                H2SO4_mass:units = "molecules/cm3air" ;
                H2SO4_mass:standard_name = "mass_concentration_of_H2SO4_in_ambient_aerosol" ;
                H2SO4_mass:long_name = "mass_concentration_of__H2SO4_in_ambient_aeroso" ;

Could you check the CMIP6 data and provide H2SO4 _mass data, please?

Thank you in advance.
Makoto Deushi.

[thomasaubry]

Dear Makoto,
Could you please confirm that:
i) It's not just a confusion around the variable name? What was called H2SO4_mass in CMIP6 was in unit of molecules/cm3air, so it looks like it was the number density. No matter how one calls it it's the same unit as what we now provide as number density.
ii) If you need a different variables, could you provide more details on what exactly you need? (the unit would have to be different from molecules/cm3air, otherwise I don't know what else than the number density that could be)

Thank you,

Thomas

[mdeushiMRI]

Dear Dr. Aubry,

I apologize for my lack of explanation.

Figure 1. H2SO4_mass (molecules/cm3air) in the CMIP6 data (multiple_input4MIPs_aerosolProperties_CMIP_IACETH-SAGE3lambda-3-0-0_gn_185001_201412.nc) at JAN2000.

Figure 2. nd (cm-3) in the CMIP7 preliminary data (nd_input4MIPs_aerosolProperties_CMIP_UOEXETER-CMIP-1-3-1_gnz_175001-202312.nc) at JAN2000.

As shown in the above figures, "nd" (cm-3) is different from "H2SO4_mass" (molecules/cm3air) by ~11 orders of magnitude.
Therefore, I think "nd" (cm-3) is a different variable from "H2SO4_mass" in the CMIP6 data (molecules/cm3air).

I think "H2SO4_mass" is defined as the number density of H2SO4 molecules contained in the stratospheric aerosols. It is the number density of NOT "aerosol particles" BUT "H2SO4 molecules".
"H2SO4_mass" (molecules/cm3air) is readily converted from the mass concentrations (in kg/cm3) of H2SO4 in the stratospheric aerosols.

Please confirm the definition noted in the CMIP6 netcdf file repeated below.

       float H2SO4_mass(latitude, altitude, time) ;
                H2SO4_mass:units = "molecules/cm3air" ;
                H2SO4_mass:standard_name = "mass_concentration_of_H2SO4_in_ambient_aerosol" ;
                H2SO4_mass:long_name = "mass_concentration_of__H2SO4_in_ambient_aerosol" ;

The "H2SO4_mass" field is necessary for our model to calculate the heterogeneous chemical reactions on the supercooled ternary solution droplets (STS, type 1b PSC).

Thank you in advance.
Makoto Deushi.

[thomasaubry]

Hello Makoto,

Apologies for the late reply. So my understanding is that you need a variable giving you the number of molecules of H2SO4 per cm3 of air, correct? If so not a problem, I will include it in the CMIP7 dataset to be released next week. I might email you test files thursday or friday to ensure it has what you are after before releasing on ESGF.

All the best,

Thomas

[mdeushiMRI]

Dear Thomas,

Thank you for your reply.

Exactly. We need the number of molecules of H2SO4 in aerosol per cm3 of air.

Best regards,
Makoto

[thomasaubry]

Perfect, thanks Makoto! I got a bit delayed but will add it over the weekend, it might then take a couple days to land on ESGF; just email me if you want a copy before.

[vnaik60]

Dear Tom (@thomasaubry)

At GFDL, we have been looking at the UOEXETER-CMIP-1-3-0 version (analysis of v1-3-1 is pending) and a question about the pre-satellite data has arisen. Attached is a plot from my colleague Fabien Paulot. My colleagues Jing Feng and David Paynter will provide more analysis later this week, but in the meantime here is an initial concern -

"I ran a 20-year running mean on the whole time series. As you can see below, it's striking that all Post Pinatubo '20-year' means are among the lowest in the entire record.
In fact, the post Pinatubo period accounts for 11 of the lowest 12 20-year period in the entire record.
Is it indeed the case that we live in an extraordinarily quiet period, or could there be some issue with this timeseries?
Thanks. Fabien"

Thanks!
Best,
Vaishali

[thomasaubry]

Thanks Vaishali! I include some plot below for SAOD anomaly (not including the not volcanic background in our dataset) and for the global mean SAOD. To help explain my reply, I've also made SAOD timeseries for the pre-satellite era only using eruptions injecting <3 Tg SO2 and <10 Tg SO2 (for reference the largest eruption of the satelitte era post Pinatubo is 2Tg SO2). Two main points:

1. If you consider all eruptions, the satellite-era period is indeed low in terms of volcanic forcing. This is not a surprise, Pinatubo is not a particularly big eruption and the 19th century was very active volcanically. I don't think there is any controversy here and our ice-core datasets should be relatively reliable in capturing past large eruptions (and we don't have other good sources of information). So yes in terms of large eruption, the satellite era is quiet.

2. If you are specifically concerned about the post-pinatubo period and the representation of small-magnitude eruptions in the dataset, have a look at the panels for eruptions <3Tg SO2 on first figure below. Most eruptions in our datasets have masses directly constrained from ice-core, but some small-magnitude eruptions are not recorded in ice-core and we had to give them a default mass depending on their volcanic explosivity index (VEI). These default masses were chosen so that the mean SAOD anomaly associated with <3Tg SO2 eruption for 1850-1979 is equal to the mean SAOD anomaly for 1998-2023, our best estimate for the forcing associated with small eruptions. So in terms of mean forcing for small eruptions, we actually designed the dataset specifically to have a zero bias. One caveat is that the geological (VEI) record is also affected by under-recording, so that in practice we lack small-magnitude eruptions, and their individual forcing on average is likely a bit too big. You can see that on the second figure below, especially in the bottom panel where you will notice that the purple curve (eruption <3 Tg SO2 pre-satellite) has typical peaks a bit higher than observed post-pinatubo (1998-2023), but also a lot of year with zero anomaly which pretty much never happens over the satellite era. You can also see the yearly SAOD distribution on the bottom right panel of figure 2 in our rough documentation (https://docs.google.com/document/d/1blX5kv0We1BteqWzMKs0OuhazAcAonay) which shows the same point. But we are talking about very small forcing signal, and again we think we did our best to avoid any bias in terms of time-average forcing for small eruptions.

Does this make sense?

[fengzydy]

Dear Tom @thomasaubry,

As Vaishali @vnaik60 mentioned earlier, David Paynter and myself at GFDL are converting the optical property data to formats required by GFDL models. When checking through SAOD data produced from different data versions (attached below), we find that the SAOD increases quite substantially from CMIP6 to v1-3-0, and we are concerned about the climate impact it may induce. For example, at 550nm (converted to GFDL band 500 to 600 nm), the average SAOD from 1850 to 1960 of cmip6plus v1-3-0 is 0.0163, it is almost twice as large as that in v1-1-3 (0.0084). Given this large increase, that would most likely double the impact on temperature and ocean heat uptake of volcanoes compared to CMIP6 (or version 1-1-3), do you have high scientific confidence that the absolute values of the volcanoes in 1-3-0 can be justified??

Thanks,
Jing

[thomasaubry]

Hello Jing,

Thanks a lot for further feedback, this is very useful!

1. I don't understand why we should remove these values or what inconsistency you refer to exactly. Both the satellite and pre-satellite era have small extinction values above these altitudes. The contribution of extinction above >=37.5 km to total SAOD is higher pre-satellite era (see bottom panel of plot below), likely because of a bias in the model we use, but GloSSAC (used for satellite era) has a lot of NaN values at such altitudes which don't necessarilly reflect a zero SAOD. In any case, should your model be unable to handle ext values above 37.5km, I suggest simply zeroing the values in our dataset. The fractional contribution of ext above 37.5km to the total SAOD is ~0.05% so I'd hope this does not affect historical simulations in any way. Such pre-processing step should simply be documented in the paper documenting your historical runs. Does this make sense/do you have any concern with the suggestion?

2. I'm not sure you will consider this as a scientific reason, but the reason is that we do not use the same aerosol model as CMIP6 :) So the whole dataset is different, really. This also affects the satellite era as all we have from GloSSAC is extinction at 525 and 1020nm. We use a reduced-complexity aerosol model (EVA_H) and, in comparison, the model (AER2D)/aerosol optical property codes used in CMIP6 are more complex and much more based on (micro)physical/chemical process. The rationale for the choice of our model is explained in the preliminary documentation document. In short, EVA_H is computationally inexpensive facilitating production of emission-derived dataset pre-satellite, its use makes us more consistent with PMIP/VolMIP, and there are major uncertainties is stratospheric aerosol models (see figure below from Clyne et al., https://doi.org/10.5194/acp-21-3317-2021, where all models are run with same SO2 injection). There is thus no immediate justification that using a complex aerosol model would produce a better forcing dataset than using with a reduced-complexity model calibrated against satellite era. Consequently, we chose an approach that maximized consistency with other MIPs, consistency between our emission and aerosol optical property products, and minimized computational cost. Hope this makes sense! The change in asymettry factor is definitly interesting and something we should comment on in documentation papers, thanks for pointing it out!

[fengzydy]

Hi Tom,

Thank you for the detailed explanation—it really helps me better understand these updates. We are currently running some experiments that reveal interesting differences in radiation and temperature fields. It’s still a work in progress, but we are trying to determine whether these changes indicate potential improvements.

Regarding extinction, I was initially concerned about the impact of small extinction values when artificially combined with thick model layers near the top of the atmosphere. Based on your response, I implemented a layer-thickness-weighted extinction as the model input, ensuring that the small values near the TOA would have minor effect.

[thomasaubry]

Thanks Jing! This sounds sensible. I'm sorry that the updates are not documented yet in a peer review publication, the CMIP timeline and the amount of work needed to make the dataset made this impossible. But writting ongoing and manuscript are hopefully out this summer :)
Looking forward to seeing simulations!

[fengzydy]

Hi Tom,

Thanks! This publication will be an important reference for CMIP7, I’m looking forward to seeing it soon.

Regarding the small extinction values in the upper atmosphere: while the small extinction values themselves are fine, we noticed that the single scattering albedo (SSA) and asymmetry factor at the same layer are significantly lower than in the rest of the stratosphere. This may lead to anomalous radiative absorption, as absorption is proportional to ext * (1 - SSA). Along those lines, we’d also like to suggest that for SSA and asymmetry factor, please replacing the small numbers at high altitude, as well as the NaNs, with nearest layer in the same band.

If SSA is not filled with valid values, users may need to take extra caution when vertically interpolating to the height/pressure layers of their models. Initially, we overlooked this detail and fill NaNs with zero, and the interpolation of SSA resulted in small SSA near the boundaries, causing quite noticeable effects on stratospheric and tropopause temperatures. In case this update could not be made, we’d be happy to provide an example interpolation script that other models could reference to avoid potential issues.

[thomasaubry]

Hello Jing, thanks for bringing this to my attention and apologies for the late reply. I'll quickly discuss the NaN value problem with task team leadership and my colleague coordinating curation of the satellite-era dataset, and post an update here.

[mo-garethsjones]

Hello Dr Aubry (@thomasaubry)
I have been looking at the v2 dataset as part of reviewing the ancillaries being created for HadGEM3/UKESM, and have been refering to the really helpful prelimiary documentation. I would like to understand why there are "NaN" values in the data. How are these to be interpreted (I see this has been discussed in a comment above)? Normally missing data values are used if the data is to be ignored.
Thanks
Gareth

[mo-garethsjones]

Hello Dr Aubry (@thomasaubry),
In the preliminary documentation it is stated that there is a contribution from "anthropogenic aerosol emissions" in the stratospheric aerosol properties. I am involved in DAMIP, and we need natural only forcing datasets for the hist-nat experiments. Is there a dataset available that has anthropogenic contributions removed?
Thanks
Gareth

[thomasaubry]

Hello @mo-garethsjones ,

I could produce a dataset that only includes the background stratospheric aerosol (i.e. without explosive volcanic eruptions). This background would still include a natural component. I'd be happy to discuss this further but I won't be able to prioritize this until picontrol, historical and scenario forcings are produced for good.

[mo-garethsjones]

Hello Dr Aubry (@thomasaubry),
if the pre-satellite era data is only constructed from emissions from a model using ice-core and geological sources, has there been any validation/comparison against available (albeit limited) optical measurements of the stratosphere pre 1970s?
Thanks
Gareth

[thomasaubry]

Hello Gareth,
To date, we've only adjusted the 20th century deposition flux in ice cores to improve the match with the pyrheliometer measurements used in CMIP6. I will document this in details in the upcoming paper. Further validation will be the focus of future work, noting that data before the 1970 is challenging to use given sparse spatial and temporal coverage, indirect and uncertain link to e.g. SAOD, and detection threshold questions. Even the 1980s are far from perfect in the satellite era! In addition to best effort on validation, future work will also include uncertainty products.
Apologies for the unsatisfying reply!

Thomas

[mo-garethsjones]

Hi Thomas, thanks for the response. After closer inspection, I am now concerned about the spatial coverage of the dataset following the Mt Agung 1963 eruption and inconsistency with pre-satellite observations. I will put my thoughts in a comment below.
Ta
Gareth

[mdeushiMRI]

Dear Dr Aubry @thomasaubry ,

As has been discussed above, we too would like to know how "NaN" values in the data are to be filled.
It would be most appreciated if you could provide the data filled with valid values.

Thanks
Makoto

[thomasaubry]

Hello @mo-garethsjones and @mdeushiMRI,
Apologies for the late reply, I'm not ignoring you! We are currently quality controlling a new version and CMIP leadership is discussing whether to release it (because of dependencies with ozone dataset). I expect a new version without these NaNs will be released but will update the github and answer questions once this has been confirmed and the files are on their way to ESGF.

[jhlee00]

Dear Dr. Aubry,

We are also currently facing issues in our KACEv2 due to NaN values in the dataset. We're filling them with zeros for now but we’re wondering if it would be better to treat them explicitly as missing values instead. Any advice on best practices would be greatly appreciated. Also, could you let us know when the filled data(new version) will be available? We're waiting for data update to run our model.

All the best,
Jaehee

[thomasaubry]

Dear @jhlee00 et al.,
The new version was finalized last week. It should be on ESGF this week (bank holiday+ESGF server changes delayed things a bit) and once it is the CMIP office will release a communication on volcanic and ozone forcing.
Thank you for your patience with this issue and apologies that the NaN problem was not detected/corrected earlier.

Thomas

[mo-garethsjones]

Hello Thomas (@thomasaubry),
I am concerned about the stratospheric aerosol dataset's spatial distribution following the 1963 Mt Agung eruption. The CMIP7 v2.0.0 daatset distribution is symmetric across the hemispheres, while the CMIP6 dataset had most of the aerosol in the southern hemisphere. The new dataset has doubled the optical depth in the Northern Hemisphere, and reduced my more than a third the optical depth in the Southern Hemisphere.

Figure of optical depth deduced from source datasets extinction diagnostics, associated with bands overlapping 550nm, monthly means for the globe, northern hemisphere and southern hemisphere.

This view is consistent with Figure 4 in https://docs.google.com/document/d/1blX5kv0We1BteqWzMKs0OuhazAcAonay

That most of the aerosol was distributed in the southern hemisphere is supported by a number of different observations and measurements ; stellar extinction, solar extinction, in-situ measurements (even some U2 flights!), and is consistent with lunar eclipse observations. Stothers, 2001, Major optical depth perturbations to the stratosphere from volcanic eruptions: Stellar extinction period, 1961-1978, JGR, has a nice summary of the evidence but there are plenty of studies from the 1960s onwards that the stratospheric aerosol impacts on radiation transmission was very much greater in the southern hemisphere. I have not seen any evidence to challenge that.

I am worried that using a dataset that potentially has too much aerosol in the stratosphere in the years following 1963 would cause unintended effects in model simulations. I had a look at the CMIP6 hist-nat experiment and examined Northern hemisphere, Southern Hemisphere and the difference between the two for TAS (1.5m air temperature). The below figure shows that the response following Agung's eruption is slightly greater in the southern hemisphere relative to the northern hemisphere, in contrast to the other major volcanic eruptions.

Figure of CMIP6 hist-nat, annual mean temperatures for Northern Hemisphere, Southern Hemisphere and the hemispheric difference.

If the more symmetric distribution is used, it will cause the simulations to have very similar responses as to the other major eruptions, with the Northern hemisphere cooling more than the southern hemisphere. Not only will this mean the spatial temperature response will be different, there will be knock on impacts on the measured impacts in models on the AMO, AMOC and others. e.g., Menary et al, 2020, Aerosol‐Forced AMOC Changes in CMIP6 Historical Simulations, GRL. Agung erupted around the same time as anthropogenic tropospheric aerosols were becoming very important to the climate response. An extra negative forcing in the northern hemisphere at this key time will add to the aerosol cooling influence in models.

Is there anything that can be done to correct for this potential bias in the aerosol distribution following Agung? I know uncertainties are great before the satellite era, but the Mt Agung eruption was a focus of a lot of scientific enquiry - the most of any large volcanic eruption up to that time. It would be a shame to not use that information to provide the best available dataset for use in models.
Thanks
Gareth

[thomasaubry]

As mentioned in the CMIP IPO communication last week, this problem has been adresses in v2.2.1. I've updated the Google doc documentation but need to update the GitHub description, thanks for the prompt. v2.0.0 will not be the recommended version. v2.2.1 is on its way to ESGF.

[jhlee00]

Dear Dr. Aubry(@thomasaubry),

In the current version(v2.2.1.) of the data, there are missing values at lower altitudes which were not present in the corresponding CMIP6 dataset. Would it be possible to fill these missing values as was done in CMIP6? Is this the final version of the dataset or will there be further updates?

Thanks
Jaehee

[thomasaubry]

Dear Jaehee,
Any remaining missing values should be for altitudes deemed tropospheric in our dataset, and as such are proper missing values. My suggestion would be to zero any missing value if you need to fill in the troposphere. Does this sound like an issue and would it be helpful if I clarified this in the preliminary documentation?
Although to my knowledge the current version (2.2.1) is not yet on ESGF owing to server changes, it's been pre-uploaded and as per CMIP IPO communication, no further update will be issued for CMIP7 AFT. So this is the final version for AFT indeed, only retractions are possible now. Tagging @durack1 @vnaik60 in case you have anything to add.

[mo-garethsjones]

Hi Thomas, I have been able to look at the v2.2.1 dataset now available in the ESGF. For completeness I include an update to the figure I put here a few weeks ago, to show the distribution of stratospheric aerosols following the 1963 Mt Agung eruption is now asymmetric across the hemispheres, with more aerosol in the Southern hemisphere than the Northern hemisphere.

Thanks for updating the dataset.

[vnaik60]

@thomasaubry you summarized the status of CMIP7 volcanic dataset well.

[durack1]

@thomasaubry apologies a late starter here.

Direct links to the datasets are below

source_id	version	ESGF direct link
UOEXETER-CMIP-2-2-1	2.2.1 (latest)	https://esgf-node.ornl.gov/search?project=input4MIPs
UOEXETER-CMIP (all versions)	1.1.3 through 2.2.1	https://esgf-node.ornl.gov/search?project=input4MIPs&versionType=all

[mingxuanwupnnl]

Dear Tom (@thomasaubry ) and all (@vnaik60 @lkemmons),

I'm Mingxuan Wu, an Earth Scientist working at PNNL. We (@hwangacme) have recently used CMIP7 volcanic SO2 emission in our US DOE E3SM earth system model. Our preliminary results show that the stratospheric sulfate burden and stratospheric AOD significantly increase by 50 to 80%. For example, SO4 stratospheric AOD and burden increase from 0.098 and 4.97 TgS to 0.157 and 8.07 TgS in 1991-12. We are concerned that this would significantly reduce surface temperature in model. We are wondering people from other modeling centers have seen this (especially for models using volcanic SO2 emission)?

[thomasaubry]

Hello Mingxuan,

So cool to see you use the SO2 inventory and running interactive :) Interactive models have a very wide range of predicted SAOD/forcing for given SO2 injections, see Clyne et al (2021, https://acp.copernicus.org/articles/21/3317/2021/). So focusing on Pinatubo, the questions are:

1. I assume you did but did you double check units? The sulfur masses we provide is Tg SO2, best to be sure you didn't think it was Tg S.
2. Were you running interactive before and if so did you have a similar injected mass for Pinatubo (15 Tg SO2)? If the answer to both of these questions is yes, unless you did significant changes to the interactive stratospheric aerosol scheme, then there is definitly an issue either with our inventory or with your implementation of it.
3. If you were running interactive before but using a different mass for Pinatubo before no suprise. If you were not running interactive before, I suspect your model is at the high end of models in terms of predicted SAOD for given mass of SO2. E.g. UKESM and CESM-WACCM need to use ~10 Tg SO2 to match Pinatubo. 15 Tg SO2 is what we deemed the best estimate but there are questions about whether these models overpredict SAOD, or whether missing processes (e.g. ash scavenging) mean the effective SO2 mass should be smaller.

I hope thinking about these questions for Pinatubo help with the wider dataset. If you suspect 3 is the "issue", you have essentially two options:
Option 1 = you could modify our emission inventory for certain eruptions to ensure that the forcing you get is something you feel comfortable with (e.g. something closer to our prescribed aerosol optical property dataset with which most modelling group will run). We strongly recommend that you document clearly all deviations from our emission inventory to facilitate inter-model comparison
Option 2 = you could use our emission inventory with no modification, which means that you might end up with a very different forcing from models running with prescribed aerosol optical properties, BUT that your simulations will be comparable to simulations from other models which run emission driven with interactive stratospheric aerosols.

For what it's worth, I think option 1 might be reassuring, but the reality is that with an interactive model you will always end up with big differences compared to our prescribed forcing. Option 2 will greatly facilitate doing good science by having simulations directly comparable between modelling group with interactive stratospheric aerosols.

Let me know what you think/any further question.

Thomas

[vnaik60]

@mingxuanwupnnl, thanks for sharing the E3SM experience with the CMIP7 volcanic SO2 emissions. It would be great if other modeling centers with interactive volcanic emissions share their experiences (ping tilmes)

Thanks @thomasaubry for providing insights and options for modelers. Our CMIP7 simulations will still use the prescribed SAOD, but for our model configuration with interactive volcanic emissions we use your option 2 as it allows us to better understand the model sensitivities keeping the forcing fixed.

[thomasaubry]

I think the approach will be similar with UKESM, i.e. I think CMIP7 sims will use prescribed aerosol optical properties, but we will run an ensemble of interactive sims with option 2.

@mingxuanwupnnl please do update us on whether you think there is a major issue with the emission inventory, or if you think it's just likely that your model is one of the models that predicts relatively high SAOD for given stratospheric SO2 emissions. I strongly suspect the latter but it will be good to know.

[mingxuanwupnnl]

@thomasaubry @vnaik60 @durack1 Thank you all for your valuable input! We (@hwangacme @keziming) just had a discussion on this and confirmed the answers to your questions. I think our E3SM model might be sensitive to the changes in interactive SO2 emission.

1. Yes. The units should be correct. We are using Tg SO2.
2. Yes. We were running interactive before. We use a volcanic emission of 10 Tg SO2 in CMIP6 and 15 Tg SO2 in CMIP7 for Pinatubo (1991-06-15).
3. When using CMIP6 SO2 emission, our model matches the sulfate burden observations well. The peak values of SAOD in our E3SM model are a bit lower than prescribed CMIP6 SAOD (when using CMIP6 interactive SO2 emission) but are a bit higher than prescribed CMIP7 SAOD (when using CMIP7 interactive SO2 emission)

@thomasaubry for your EVA model, I wonder what the configuration of the gas chemistry it is. Does it have interactive stratopsheric ozone and other oxidants (e.g., OH).

[thomasaubry]

Thanks for confirming @mingxuanwupnnl ! No problem with the SO2 inventory then. The 15 Tg SO2 value from CMIP7 is simply the one recorded in the satellite era SO2 emission inventory we use (MSVOLSO2L4, curated by Simon Carn/NASA), and it also falls well within the range of estimates for Pinatubo (e.g. Guo et al., 2004). The discrepency between some model simulations when using SO2 masses >15 Tg and observations in terms of sulfate burden or SAOD is well known (e.g. Timmreck et al., 2018; Qaglia et al., 2023). Again this raises the question of missing processes in the models (ice/SO2 scavenging) often invoked to justify injecting 10Tg SO2 as an "effective" mass, or of other model biases (some models predicts much lower SAOD than others for same SO2 mass, see Clyne et al., 2021, Qaglia et al., 2023). Our role as dataset providers is not to provide the SO2 mass that will make most models match aerosol burden/SAOD observations, but to provide the SO2 mass that is consistent with SO2 observations (noting of course observations for SO2 and SAOD are uncertain themselves). This brings you back to option 1 vs option 2 in my September 11 post :)

An extensive dataset documentation is coming soon in GMD but as pointed in our initial short documentation, EVA_H is a reduced complexity volcanic aerosol model, not an interactive stratospheric aerosol model. The version used in CMIP7 was calibrated using the CMIP7 1979-2023 SO2 inventory as input to match the CMIP7 stratospheric aerosol optical property dataset. The original model EVA was widely used to generate forcings for some CMIP6 MIPs (e.g. PMIP, VolMIP) and other applications because of the large discrepencies among stratospheric aerosol models discussed in the paragraph above and references therein.

I hope this clarifies!

[MartinDix]

Our model uses AOD. https://docs.google.com/document/d/1blX5kv0We1BteqWzMKs0OuhazAcAonay gives the 1850-2021 mean as 0.0135 and my calculation reproduces this. However there is a small inconsistency in the climatology files where I calculate a mean of 0.0137.

This seems to be due to a change in the treatment of the tropopause (level below which data is missing) after 1979. Before this the tropopause is constant in time and afterwards it varies seasonally. I can match the climatology values by calculating the mean over non-missing values Taking missing values as zero gives the global value of 0.0135 which seems more correct.

[thomasaubry]

Dear Martin,

Sincere apologies for the late reply! I will double check I find the same issues but I recently found a similar problem with my Scenario files so I suspect that's the case indeed. I'll see if we need to suggest an implementation recommendation (likely following yours) here, in the google doc and in the revisions of the documentation paper. But PI and historical datasets will not change for sure at this stage, and the inconsistency is luckily quite small. Thanks for spotting this!

[ilaria-quaglia]

Hi @thomasaubry,
I’m currently preparing the emissions file for the scenario simulations using the 1950–2021 emissions dataset. In your preprint (https://doi.org/10.5194/egusphere-2025-4990), you mention that for the stratospheric aerosol climatology a 9-year ramp over 2022–2030 is applied for scenario simulations to ensure a smooth transition from the historical period.
I was wondering whether a similar ramping approach should be applied to the SO₂ emissions file, or if instead I should use the same SO₂ file that is used for the pre-industrial period.

In addition, while looking at the code you provide to calculate the climatology (https://github.com/thomasaubry/CMIP7_stratforcing_v2.2.1/blob/main/makeclimatology.m), I noticed that the climatology contains 12 distinct monthly values. I have not included any seasonality and instead applied the same annual mean value to each month. Do you think this approach is still appropriate?

[thomasaubry]

Hello Ilaria,

Sorry that we haven't provided PI and Scenario emission files! Thanks for your questions and I'm always happy to be contacted to support creation of PI/Scenario emission files.

My recommendation would be to be as consistent as possible with the aerosol optical property dataset, so indeed have a linear ramp between Dec 2021 and the PI climatology starting Jan 2031. For the same reason, I'd also recommend including the seasonality in the PI file, but if you started running PI simulations with seasonally-invariant emissions, I would likely not worry too much about it (just document it in the paper documenting sims). The most important seasonal component is likely the one in the non-volcanic background, which your model will simulate interactively.

Last, did you mean 1850-2021 instead of 1950-2021?

Thanks!

Thomas

[ilaria-quaglia]

Thank you! yes, sorry, I meant 1850 :)

[DirkOlivie]

Hi @thomasaubry

we are in the process of preparing the volcanic forcing data (ext, ssa, asy) for the NorESM3 model based on the CMIP7 files from input4MIPs (version UOEXETER-CMIP-2-2-1).

My impression is that the optical properties in the climatology files (ext/ssa/asy_input4MIPs_aerosolProperties_CMIP_UOEXETER-CMIP-2-2-1_gnz_185001-202112-clim.nc), are calculated as a standard average based on the time series file. I however would have expected that a weighted average should have been used to calculate the time average of ssa (where one uses ext as weighting) and the time average of asy (where one uses ext*ssa as weighting).

II think that in the volcanic forcing data made available for CMIP6, such a weighted average had been used to generate the mean climatology of the optical parameters.

Is there any reason why ssa and asy are standard averages and not weighted averages in the CMIP7 climatology? Is it best that modelling centers use the climatology provided in UOEXETER-CMIP-2-2-1, or is better to use a weighted averaged climatology (which the modelling centers calculate themselves)?

Best regards,
Dirk Olivié

[DirkOlivie]

Hi @thomasaubry, @vnaik60, @thomasaubry,

I have now run 4 simulations with NorESM testing different volcanic forcing climatologies. The simulations are :
(S1) simulation without any volcanic forcing,
(S2) simulation with CMIP6 volcanic forcing climatology (weighted average over 1850-2014 period),
(S3) simulation with CMIP7 volcanic forcing climatology provided by CMIP7 (standard average over 1850-2021), and
(S4) simulation with CMIP7 volcanic forcing climatology derived from the CMIP7 time series (weighted average over 1850-2021).

The simulations have a length of 5 years, and the horizontal wind is nudged to ERA5 (to increase the signal-to-noise ratio when comparing radiative fluxes in the above 4 simulations).

In the table below are shown global mean radiative fluxes [W/m2] averaged over the 5 years. The first three data columns show differences with respect to the reference simulation (S1, without any volcanic aerosol); the last data column shows the difference between simulations (S3, CMIP7) and (S4, CMIP7-weighted). The lowest line in the table shows the AOD (in the 425-625 nm band). Positive TOA fluxes are downward.

	CMIP6	CMIP7	CMIP7-weighted	CMIP7 minus CMIP7-weighted
	(S2)-(S1)	(S3)-(S1)	(S4)-(S1)	(S3) - (S4)
SW+LW at TOA	-0.207 (+/-0.004)	-0.280 (+/-0.008)	-0.231 (+/-0.007)	-0.049 (+/-0.007)
SW+LW at TOA (clear-sky)	-0.311 (+/-0.004)	-0.429 (+/-0.005)	-0.354 (+/-0.007)	-0.075 (+/-0.004)
SW surface downwelling	-0.211 (+/-0.023)	-0.301 (+/-0.021)	-0.269 (+/-0.029)	-0.033 (+/-0.013)
SW surface downwelling (clear-sky)	-0.411 (+/-0.008)	-0.591 (+/-0.011)	-0.525 (+/-0.014)	-0.066 (+/-0.005)
AOD volcanic (425-625 nm)	0.010	0.012	0.012	0.0

So the use of (S3) instead of (S4) leads to a difference in the net TOA flux of around -0.05 W/m2 (cooling). At the surface the SW downwelling radiation is weaker (-0.033 W/m2, cooling) when using (S3) compared to (S4).

[durack1]

@DirkOlivie thanks so much for pushing through these simulations, it's wonderful to have these direct comparisons to ponder - and the magnitude is eye-popping.

Just so I can catch up, are the tabulated values about the area-weighted global mean values? If yes, I can imagine there is a rich spatial pattern that also emerges, giving rise to the large differences (~0.1) when comparing clear-sky. Also, where does NorESM3 sit in the range of CMIPx aerosol sensitivities? I know that CESM and E3SM were sensitive, so attempting to do a back-of-the-envelope calculation, would this ~0.05 W m-2 be small/average/large compared to other model sensitivities if you had to guess?

[DirkOlivie]

Hi @durack1, thanks for your reaction.

The tabulated values are indeed area-weighted global mean values. To see the structure of the differences, I have now made some zonal-averaged plots and some 2D maps (see below). The test simulations have been done with NorESM2.3 (in CMIP6 we used NorESM2.0, in CMIP7 we aim to use NorESM3.0) - NorESM2.3 is an intermediate version. The NorESM2.3 model has indeed a strong negative aerosol effect (mainly from aerosol-cloud interactions), but I assume that it does not play a large role for the volcanic aerosol forcing (which goes mainly through aerosol radiation interaction). So the numbers obtained here might possibly be representative also for other models.

The plots below show the zonal mean distribution of the fluxes in the table above.

The plots below show maps of the clearsky TOA SW+LW flux (left), and of the clearsky SW downwelling flux (right). Shown are differences for CMIP6 w.r.t. no forcing (upper row), CMIP7 w.r.t. no volcanic forcing (middle row) and CMIP7-weighted w.r.t. no volcanic forcing (bottom row).

[znichollscr]

Thanks @DirkOlivie! I've added a pointer to this conversation into the docs so that others can more easily find it, see #440

[thomasaubry]

thanks a lot for running the sims and making the analysis available here @DirkOlivie. Not a massive difference but not small either, definitly wish we had caught that before.