Tuning Leung_2023 dust cycle for clm6_0_cam7.0 and namelist changes

[dmleung]

In issue #1423 we asked for dust tuning in CAM for CESM3 because dust lifetime is generally too low and it requires really high dust emissions to sustain a reasonable dust AOD. This PR fixes issue #1423.

This PR adds in two fixes relevant to coarse-mode aerosol asphericity on aerosol lifetime and aerosol optics. The third fix is related to dust lifetime by tuning the dust size distribution from dust emission.

1. Observational and theoretical studies showed that aspherical coarse-mode dust enhances dust optical thickness and mass extinction efficiency (MEE) by 20–60 % (Fig. 1d of Jasper Kok et al., 2017). We add in a parameter dustaspherical_opts that scales up dust AOD by an extra 30 %.
2. Observational and lab studies showed that aspherical coarse-mode dust has reduced gravitational settling velocity by 20 %. We add in a parameter asphericaldust_drydep that scales down the gravitational settling velocity of coarse-mode aerosols by 20 %.
3. We change the mass fraction of emitted dust size distribution dust_emis_sclfctr, by increasing the mass fraction of the emitted accumulation-mode dust from 1% to 4 %, and reducing the mass fraction of the coarse-mode dust from 98.9 % to 96 %.

We eventually add changes to dust_emis_fact to make sure Leung_2023 global dust AOD stays around 0.029 for the 2000s in the F cases.
With all these changes, we allow global mean dust AOD to be around 0.029 but substantially reduce dust emission down from 5.3 Pg/yr to 2.5 Tg/yr. This means CAM now has a more reasonable d(dustAOD)/d(dustemis) sensitivity.
This dust emission budget agrees much better with the CMIP6 multimodel agreement of 1.5–3 Pg/yr.

Note:

1. The tuning parameter for dust optics (dustaspherical_opts) is a quick fix from Longlei Li @L3atm. In the long-term we may introduce a modification to the Mie theory that explicitly calculates the mass extinction efficiency (MEE) of aspherical dust and other aerosols.
2. For the gravitational settling and dust size distribution, the adjustments are from empirical results from papers and there are no equations. These parameters will likely stay around for tuning purposes.
3. All the above changes directly apply to both Zender_2003 and Leung_2023 dust.

Contributors: @dmleung @L3atm @tilmes @jfkok

@fvitt @dlawrenncar @ekluzek

[fvitt]

@dmleung I believe the code mods for the dust optics should be done in refractive_aerosol_optics_mod.F90, rather than aerosol_optics_cam.F90.

[dmleung]

@fvitt maybe we can talk about this if you have an idea how to do this in refractive_aerosol_optics_mod.F90, thanks

[dmleung]

@fvitt @tilmes
After some study on the refactored optics modules, it seems like the function refractive_index_sw uses the refractive indices with the volume averaging method to sum up all refractive indices across all aerosol species together:
do icol = 1, ncol
vol(icol) = specmmr(icol,ilev)/specdens
crefin(icol) = crefin(icol) + vol(icol)*specrefindex(iwav)
end do
This seems to be the only occasion where we can do something with a single species. Later on, crefin is used to calculate extinction cext and pext in refractive_aerosol_optics_mod inside sw_props. However, since refractive index does not scale linearly with extinction efficiency cext from the Mie theory, I cannot just scale up crefin by a certain amount here.

Alternatively, maybe we can change the look up table for dust directly. I am still studying the code on the look up table (terms like extpsw).
See the plot below that illustrates the difference between spherical optics using the Mie theory (blue) and the aspherical optics (orange). The discrepancy between blue and orange in the coarse mode range is why we want to scale up extinction for coarse dust by 30 %. Can we put these values into the look up table? (This curve only works for dust, we can change something to make it work on seasalt too.)

@tilmes may also have an opinion on how to change things using this plot.

[dmleung]

@fvitt @tilmes
I have done some code cleanup. In the current setting, we set optics and deposition changes to apply to MAM and not CARMA/BAM. I did this by using a flag called modal_active. (We can add carma_active too if we want to apply the changes to CARMA.) For modal_active=.True., I mainly added ntot_amode that indicates how many modes there are for different MAM versions. Then, we use n_coarse_dust to indicate which mode is the coarse dust mode to apply dust asphericity-induced optical and aerodynamic changes.

If we want to make the changes effective for CARMA, we will need to figure out how to indicates bins with D > 1 μm.
I also did code cleanup, and I think Francis can start reviewing the code (and, if desired, make necessary changes to make CARMA work too). Thanks.

[dmleung]

@fvitt @tilmes The CAM dust tuning is done and ready for review. The science changes are all done. I merged the cam tag to the latest 6_4_128.
I also made changes in the namelist_defaults_cam.xml for dust_emis_fact. Note that we support dust_emis_fact values for
CAM7 with Leung_2023 and CAM6 with Zender_2003. We don't have dust_emis_fact values for neither CAM6 with Leung_2023 nor CAM7 with Zender_2003.
After discussions with @dlawrenncar, @wwieder, and @ekluzek, we decided for now we do not prescribe separate dust_emis_fact values for cases with CTSM in the BGC mode (e.g., the B cases).

The following plot shows dust from Zender_2003 in CLM5/CAM6 (left) and Leung_2023 in CLM6/CAM7 (right) in the 2000s for the FHIST and the FHIST_LTso cases, with CTSM using the SP mode. This should be how dust in CESM3.0 look like.
We successfully scaled Leung_2023 dust emission down from > 5000 Tg/yr (which caused the issue of super high surface dust PM10) to < 3000 Tg/yr while maintaining roughly the same levels of dust extinction. This is ready for @cecilehannay and @adamrher to do their tuning runs.

In the future we will try to continue working on improving the CAM aerosol lifetime.

[dmleung]

After testing the dust_emis_fact in the FHISTC_LTso case (SP mode with climo LAI) without nudging, the following dust_emis_fact values work best.

These two will be in namelist_defaults_cam.xml:
Leung_2023 in CLM6/CAM7 physics: dust_emis_fact = 4.0D0 (2.3D0 before tuning)
Zender_2003 in CLM5/CAM6 physics: dust_emis_fact = 1.75D0 (0.7D0 before tuning)

This increase in dust_emis_fact means that we reduce dust emissions and burden in CAM.

The following will not be supported:
Leung_2023 in CLM5/CAM6 physics: dust_emis_fact = 8.0D0
Zender_2023 in CLM6/CAM7 physics: dust_emis_fact = 1.3D0

In CESM2.2.2, online nudging (&nudging_nl) will mostly reduce global dust by ~10–20 %, so the dust_emis_fact values should be a bit smaller than the ones displayed above to compensate the dust reduction due to online nudging.

We did not make tests for CAM5 and we are not changing namelists for CAM5.

[adamrher]

Hi @dmleung @tilmes I'm wondering if we can or should do an add'l science validation comparing the new AOD to the MODIS and MERRA2 products? I think this is included in some version of ADF? The reason I ask is I'm finding it difficult to bring down RESTOM in CESM3, and I'd like to double check whether the dust AOD can or should be increased marginally to help bring down the clear-sky SW fluxes. I think it would be useful to see the comparisons in an F-case (SP) and B-case (BGC). For the B-case you can use the latest coupled case NCAR/cesm_dev#209. The seasalt_emis_scale = 1.0 in this run, and the corresponding F-case should probably also use that value.

Is that reasonable or perhaps you've already done a careful comparisons to obs?

[tilmes]

@adamrher I performed a nudged model simulation using alpha07e.ctsm5.3.075, before Danny's dust tuning, but we set the seasalt_emis_scale to 1.0.

We nudged to MERRA2 for 2016-2018 and I compared to ATom data. The run is here:
/glade/campaign/cesm/cesmdata/cseg/runs/cesm2_0/f.e30.alpha07e.ctsm5.3.075.sp.nudged.2016_2018_MTt4s.ne30.001

The comparisons are the MTt4S run (green), a different CAMchem MAM5 run with sea-salt emissions = 1 (blue) and a CARMA run (red), that we are working on now.

So, you want to look at the green line compared to the black. Here is sea-salt:

We concluded:

• seasalt is in general very good, with a little larger values at the surface, and an overestimation in the Pacific 30-60S, while there is a underestimation in the Pacific 0-30S.

We could try to increase the factor slightly, maybe to 1.2 and see how we are doing?

I can now add Danny's changes as well, is there a new tag I can use or should I use some SourceMods?

Regarding other aerosols, like sulfate:

We still have some strong underestimation in all the remote regions below 5-6km, while we have some higher values values above. We are working on new processes to include in the model to improve sulfate, but that is not for the current version.

Organic aerosols are also too low in the lower troposphere, but OK higher up.
Finally, we do still overestimate extinction mostly in the tropics over the Atlantic.

[tilmes]

@adamrher you can run the AOD comparisions with MODIS and MERRA, they are in the default ADF.. not sure how to turn on though, Justin should know.

[adamrher]

@tilmes I'll have Cecile run the ADF with the comparisons to MODIS and MERRA2. I'd like to see how we're doing in the free running runs (not nudged). I have a FHISTC_LTso run with seasalt set to 1.0 and Danny's dust mods, so I'll use that.

[adamrher]

I am nix'ing any additional science validation from holding up this PR. Sorry for the firedrill.

[dmleung]

Hi @adamrher, below is a plot on dust AOD evaluation against MODIS dust AOD (MIDAS reanalysis) in the 2000s. This run uses the latest debugged code with the compset FHISTC_LTso without nudging. I like the current dust level (with a global mean of ~0.027-0.028), which gives a better comparison against MIDAS dust AOD. But, I would be okay with pushing the dust AOD up to ~0.03 if really needed.

[adamrher]

Thanks @dmleung! I've pinged you in a google slide show to follow up on these results. While we may decide to change the dust_emis_fac proposed in this PR (4.0), I don't want that to slow down this PR so it should proceed as is with the 4.0 value.

[cacraigucar]

A couple of restart regression tests are failing with differences between results from a full run and a restart run:

See: /glade/derecho/scratch/cacraig/aux_cam_intel_20251119160153/ERP_D_Ln9_P64x2.f09_f09_mg17.QSC6.derecho_intel.cam-outfrq9s.GC.aux_cam_intel_20251119160153

In this run, the following two fields are different in a normal CAM output file:

 AODDUST   (lon,lat,time)  t_index =      1     1
          4    55296  (    27,   150,     1) (   187,    24,     1) (   151,   182,     1) (   151,   182,     1)
               34560   2.766596111081480E-28   8.196876032849101E-29 3.2E-32  1.294342741744524E-28 1.7E-08  1.294342741744524E-28
               34560   2.766596111081480E-28   8.196876032849101E-29          1.294665525454745E-28          1.294665525454745E-28
               55296  (    27,   150,     1) (   187,    24,     1)
          avg abs field values:    1.420712484107264E-28    rms diff: 2.3E-34   avg rel diff(npos):  1.7E-08
                                   1.420712506523419E-28                        avg decimal digits(ndif):  3.9 worst:  3.6
 RMS AODDUST                          2.3102E-34            NORMALIZED  1.6261E-06

 AODDUST03   (lon,lat,time)  t_index =      1     1
          4    55296  (    20,    41,     1) (   187,   101,     1) (   151,   182,     1) (   146,    81,     1)
               34560   1.111193380734803E-28   1.706131086478035E-29 2.8E-32  5.191711338743631E-29 4.9E-08  2.741491162962787E-29
               34560   1.111193380734803E-28   1.706131086478035E-29          5.194516466735632E-29          2.743428253648504E-29
               55296  (    20,    41,     1) (   187,   101,     1)
          avg abs field values:    4.443324053750224E-29    rms diff: 2.0E-34   avg rel diff(npos):  4.9E-08
                                   4.443324238612640E-29                        avg decimal digits(ndif):  3.5 worst:  3.2
 RMS AODDUST03                        1.9579E-34            NORMALIZED  4.4063E-06

On izumi (the local cluster) the following test had differences in the cpl history file. See:
/scratch/cluster/cacraig/aux_cam_gnu_20251119142346/ERP_Ln9_P24x2.f45_f45_mg37.QPWmaC6.izumi_gnu.cam-outfrq9s_mee_fluxes.GC.aux_cam_gnu_20251119142346

[cacraigucar]

A couple of restart regression tests are failing with differences between results from a full run and a restart run:

See: /glade/derecho/scratch/cacraig/aux_cam_intel_20251119160153/ERP_D_Ln9_P64x2.f09_f09_mg17.QSC6.derecho_intel.cam-outfrq9s.GC.aux_cam_intel_20251119160153

In this run, the following two fields are different in a normal CAM output file:

 AODDUST   (lon,lat,time)  t_index =      1     1
          4    55296  (    27,   150,     1) (   187,    24,     1) (   151,   182,     1) (   151,   182,     1)
               34560   2.766596111081480E-28   8.196876032849101E-29 3.2E-32  1.294342741744524E-28 1.7E-08  1.294342741744524E-28
               34560   2.766596111081480E-28   8.196876032849101E-29          1.294665525454745E-28          1.294665525454745E-28
               55296  (    27,   150,     1) (   187,    24,     1)
          avg abs field values:    1.420712484107264E-28    rms diff: 2.3E-34   avg rel diff(npos):  1.7E-08
                                   1.420712506523419E-28                        avg decimal digits(ndif):  3.9 worst:  3.6
 RMS AODDUST                          2.3102E-34            NORMALIZED  1.6261E-06

 AODDUST03   (lon,lat,time)  t_index =      1     1
          4    55296  (    20,    41,     1) (   187,   101,     1) (   151,   182,     1) (   146,    81,     1)
               34560   1.111193380734803E-28   1.706131086478035E-29 2.8E-32  5.191711338743631E-29 4.9E-08  2.741491162962787E-29
               34560   1.111193380734803E-28   1.706131086478035E-29          5.194516466735632E-29          2.743428253648504E-29
               55296  (    20,    41,     1) (   187,   101,     1)
          avg abs field values:    4.443324053750224E-29    rms diff: 2.0E-34   avg rel diff(npos):  4.9E-08
                                   4.443324238612640E-29                        avg decimal digits(ndif):  3.5 worst:  3.2
 RMS AODDUST03                        1.9579E-34            NORMALIZED  4.4063E-06

On izumi (the local cluster) the following test had differences in the cpl history file. See: /scratch/cluster/cacraig/aux_cam_gnu_20251119142346/ERP_Ln9_P24x2.f45_f45_mg37.QPWmaC6.izumi_gnu.cam-outfrq9s_mee_fluxes.GC.aux_cam_gnu_20251119142346

This looks like it could be something not being done in double precision (the differences are of the right relative size). We have an automated checker for constants, and that checker indicated all of the constants were defined with the "r8" and I looked at the namelist constants, and I didn't see any without the "D0". I looked at the other modified code, and didn't see any places where it might be doing single precision arithmetic. I am pointing this out though, as I still consider it a possibility.