Bug in jacobian background simulation stratosphere

[eastjames]

Name and Institution (Required)

Name: James East Institution: ACMG

Description of your issue or question

I ran into a bug background simulation for the jacobian runs. This simulation should have 0 emissions, but a pattern that looks like numerical instability shows up partway through the simulation. The CH4 concentrations blow up to unphysical values (~600,000 ppb)

The plots show simple average CH4 from the GEOS-Chem SpeciesConc files. The bug appears to happen in the stratosphere.

This is using lognormal errors and the feature/2d_sensitivities branch. The error is in the GEOS-Chem output and doesn't result from any IMI changes in that branch.

[eastjames]

Tagging @msulprizio

[msulprizio]

We have seen similar patterns in other nested grid simulations of extremely high concentrations in a checkerboard pattern. It often has to do with the Courant number calculation. Others have been able to get around it by increasing the transport timestep and/or by increasing the nested grid buffer cells (default is 3 cells on each edge). You may also want to check the boundary conditions and/or met fields to see if anything stands out there as triggering this issue.

See also similar issues:

• https://github.com/geoschem/geos-chem/issues/605
• https://github.com/geoschem/geos-chem/issues/2241

[eastjames]

Halving the advection time step in geoschem_config.yml solved the issue. This suggests the problem is most likely related to the CFL condition being violated, pointed out by @msulprizio and @hannahnesser in offline conversations.

Old:

#============================================================================
# Timesteps settings
#============================================================================
timesteps:
  transport_timestep_in_s: 300

New:

#============================================================================
# Timesteps settings
#============================================================================
timesteps:
  transport_timestep_in_s: 150

[nicholasbalasus]

For sake of reproducing this issue, here's a script to run GEOS-Chem with the typical time step (300 s) and the reduced time step (150 s). Plotting the maximum methane concentration in any grid box in the domain shows that the more frequent transport time step prevents this issue (as James found out). It is an imperfect comparison, but interestingly, the 300 s simulation took ~4 hours and the 150 s simulation took only ~4.75 hours.

Run script

``` #!/bin/bash for transport_timestep in 150 300; do # GC for [CH4, GEOS-FP, 0.25 x 0.3125, NA, 47L] source /n/home06/nbalasus/envs/gcclassic.rocky+gnu10.minimal.env dir="/n/holyscratch01/jacob_lab/nbalasus/for-james/transport_${transport_timestep}" rm -rf "${dir}" mkdir -p "${dir}" cd "${dir}" git clone https://github.com/geoschem/GCClassic.git cd GCClassic git checkout 14.4.1 git submodule update --init --recursive cd run run_dir="gc_run" c="9\n2\n4\n4\n2\n${dir}\n${run_dir}\nn\n" printf ${c} | ./createRunDir.sh cd "${dir}/${run_dir}/build" cmake ../CodeDir -DRUNDIR=.. make -j make install cd "${dir}/${run_dir}" # Configuration files sed -i -e "s|'Metrics',|#'Metrics',|g" \ -e "s|'CH4',|#'CH4',|g" \ -e "s|00000100 000000|00000000 010000|g" \ -e "s|time-averaged|instantaneous|g" HISTORY.rc sed -i -e "s|GEOS_0.25x0.3125|GEOS_0.25x0.3125_NA|g" \ -e "s|\$RES|\$RES.NA|g" HEMCO_Config.rc.gmao_metfields bc="/n/holylfs05/LABS/jacob_lab/imi/ch4/blended-boundary-conditions/v2024-06" sed -i -e "s|RF xy|C xy|g" \ -e "s|\$ROOT/SAMPLE_BCs/v2021-07/CH4|${bc}|g" \ -e "s|./Restarts/GEOSChem.Restart.\$YYYY\$MM\$DD_\$HH\$MNz.nc4|${bc}/GEOSChem.BoundaryConditions.\$YYYY\$MM\$DD_0000z.nc4|g" \ -e "s|SpeciesRst|SpeciesBC|g" HEMCO_Config.rc sed -i -e "s|20190101|20230101|g" \ -e "s|20190201|20230201|g" \ -e "s|-130.0, -60.0|-128.75, -61.25|g" \ -e "s|9.75, 60.0|10.5, 59.25|g" \ -e "s|transport_timestep_in_s: 300|transport_timestep_in_s: ${transport_timestep}|g" geoschem_config.yml # Run script cp runScriptSamples/operational_examples/harvard_cannon/geoschem.run . sed -i -e "s|-c 8|-c 32|g" \ -e "s|0-12:00|0-24:00|g" \ -e "s|15000|128000|g" geoschem.run sbatch geoschem.run done ```

Plotting script

``` import re import glob import pandas as pd import xarray as xr import matplotlib.pyplot as plt # For each hour, collect the maximum CH4 concentration in any grid cell dir150="/n/holyscratch01/jacob_lab/nbalasus/for-james/transport_150/gc_run/OutputDir/" dir300="/n/holyscratch01/jacob_lab/nbalasus/for-james/transport_300/gc_run/OutputDir/" files = sorted(glob.glob(dir150 + "*SpeciesConc*.nc4")) hours = [re.search(r'(\d{8}_\d{4})', f).group(0) for f in files] max_ch4_150, max_ch4_300 = [],[] for hour in hours: file = f"GEOSChem.SpeciesConc.{hour}z.nc4" with xr.open_dataset(dir150+file) as ds: max_ch4_150.append(1e9*ds["SpeciesConcVV_CH4"].values.max()) with xr.open_dataset(dir300+file) as ds: max_ch4_300.append(1e9*ds["SpeciesConcVV_CH4"].values.max()) # Plot the max concentrations for each simulation fig,axs = plt.subplots(2,1,figsize=(14,10)) axs[0].plot(pd.to_datetime(hours, format="%Y%m%d_%H%M"), max_ch4_300, label="300 s transport time step", color="#ee6677") axs[1].plot(pd.to_datetime(hours, format="%Y%m%d_%H%M"), max_ch4_150, label="150 s transport time step", color="#4477aa") axs[1].ticklabel_format(axis='y', style='sci', scilimits=(3,3)) for ax in axs: ax.legend() ax.grid(axis="y", linewidth=0.25) ax.set_ylabel("Maximum Grid Box CH$_4$ Concentration [ppb]") ```