Mass conservation issue for chemical tracers in a continuous run mode

[jianheACM]

What is wrong?

In a continuous run mode, meteorological fields are reinitialized at 00z for each 24-hour cycle and chemical fields are read from the restart files (except for the first day) and processed through a python script in the global workflow (i.e., https://github.com/NOAA-EMC/global-workflow/blob/develop/scripts/exgfs_aero_init_aerosol.py).
I add a passive tracer CO2 with initial condition of 400 ppm into the model. There is no source or sink for CO2. I let the model run in a continuous run mode as described above. The global total mass of CO2 keeps decreasing. With a 35-day continuous run, the global total mass of CO2 could decrease by 2%, which will end up with 21% in a year-long simulation. Detailed description could be found here.
@zhanglikate @rschwant @bbakernoaa

What should have happened?

The global total mass of CO2 should be conserved or at least comparable to what we see in a free run mode, which is about 0.0286% per year. See slide 7 of the document.

What machines are impacted?

All or N/A

Steps to reproduce

I started an experiment with C96 and 65-level configurations.
I add a passive tracer CO2 with initial condition of 400 ppm into the model. There is no source or sink for CO2. I let the model continuously run (meteorological fields are reinitialized at 00z for each 24-hour cycle and chemical fields are read from the restart files except for the first day and processed through a python script in the global workflow) for 20160715-20160823 and have hourly output for CO2.
I calculate global total mass of hourly CO2 and see a decreasing trend for the simulation period.

Additional information

In the global workflow, specifically in https://github.com/NOAA-EMC/global-workflow/blob/develop/scripts/exgfs_aero_init_aerosol.py, there is an adjustment factor based on the ratio of delp in the restart file and initial condition file, which is applied to chemical fields from the restart files. However, this is not a real vertical interpolation. Detailed description could be found slides 2-4 here.

Do you have a proposed solution?

I have tested a different method to process chem tracers from restart files. See slides 5-6 of the document. Results could be found on slide 7.
Despite the mass conservation issue is improved in the testing method, it is still not as good as in the free run mode. So we need a better way to improve the mass conservation issue for chemical tracers.
Alternatively, we can perform a vertical interpolation for chemical tracers in https://github.com/NOAA-EMC/global-workflow/blob/develop/scripts/exgfs_aero_init_aerosol.py instead of simply applying an adjustment factor. But this will end up with two vertical interpolations for chemical tracers, one in the global workflow and one in the model source code.

[DWesl]

You mention sophisticated functions to perform vertical interpolation while preserving mass between pressure levels. Could you elaborate on what those functions are?

[zhanglikate]

You mention sophisticated functions to perform vertical interpolation while preserving mass between pressure levels. Could you elaborate on what those functions are?

@rmontuoro is more familiar with the python script vertical interpolation script of https://github.com/NOAA-EMC/global-workflow/blob/develop/scripts/exgfs_aero_init_aerosol.py , who can help to elaborate more details about it.

[DWesl]

It looks like that script passes the processing to another:

global-workflow/scripts/exgfs_aero_init_aerosol.py

Line 52 in b405b7d

merge_script_pattern = "{ush_gfs}/merge_fv3_aerosol_tile.py"

global-workflow/scripts/exgfs_aero_init_aerosol.py

Line 286 in b405b7d

subprocess.run([merge_script, atm_file, tracer_file, core_file, ctrl_file, rest_file, tracer_list_file, temp_file], check=True)

which currently does the aforementioned rescaling here:

global-workflow/ush/merge_fv3_aerosol_tile.py

Lines 70 to 71 in b405b7d

	# read pressure layer thickness from restart file
	delp = rest_file["delp"][0, :]

global-workflow/ush/merge_fv3_aerosol_tile.py

Lines 91 to 93 in b405b7d

	dp = np.zeros(delp.shape)
	for k in range(0, dp.shape[0]):
	dp[k, :, :] = ak[k + 1] - ak[k] + psfc * (bk[k + 1] - bk[k])

global-workflow/ush/merge_fv3_aerosol_tile.py

Line 95 in b405b7d

scale_factor = delp / dp

global-workflow/ush/merge_fv3_aerosol_tile.py

Lines 112 to 113 in b405b7d

	base_file[variable_name][0, :, :] = 0.
	base_file[variable_name][1:, :, :] = scale_factor * variable[0, :, :, :]

It looks like that script rescales the CO₂ mixing ratio to keep the mass per layer constant and sets the layer not available in the restart file to zero.

Given the difference in number of levels between the restart file and the file the model reads in, would setting base_file[variable_name][0, :, :] = base_file[variable_name][1, :, :] after line 113 make sense? If the mass differences printed on lines 105-120 match the mass differences found in the experiment, the zeroed layer is probably not the problem. (EDIT: This is adding a zero on the 1 Pa end of the model, so that in itself would not cause much of a difference, though the lower-CO₂ air aloft being moved surfaceward in the model-internal-to-restart-file transition might still be a concern)

I am reminded that, if the surface pressure changes, it is possible to conserve the CO₂ mass per unit area or mole fraction/mixing ratio values but not both (that is, keeping the mole fraction values at 400ppm will result in the column CO₂ mass changing with the surface pressure).

The issue mentions a proposed solution, described in the accompanying presentation, that achieves better mass conservation than the existing method. I would like more details on this new method.

As an alternative, is there an option to have the model runs be three or five days at a stretch rather than one? That would result in fewer shocks to the system from new meteorological data, with more time as a free run, which doesn't have those problems.

EDIT: Just remembered that runs into problems. Is there an analysis-nudging option?

EDIT2: atmos_cubed_sphere has a nudging option, but global-workflow doesn't seem to have integrated it. Fortunately, it seems to prefer adjusting just temperature and winds, which should allow a simulation started with 400 ppm everywhere and no fluxes to continue forward with 400 ppm everywhere and the same total mass. Unfortunately, there doesn't appear to be any documentation on what kind of files it wants as input.

[yangfanglin]

Daniel, thank you for checking. Looks like the rescaling approach needs further consideration.

I noticed also on slide 3 of Jian He's ppt there is notation Ptop=0. Please note FV3 models cannot have a model top at zero. Ptop should be 20 Pa. If there is a zero pressure in gfs_data.tile[1-6].nc it should have never been used in any model integration. Not sure if this has played any role in interpolation and calculation of mass.

[WalterKolczynski-NOAA]

Is this even going to be a problem anymore since we are going to be using aerosol DA? The old method of continuous aerosols needing to adjust restarts from the previous forecast won't be needed, AFAIK.

[DWesl]

Is this even going to be a problem anymore since we are going to be using aerosol DA? The old method of continuous aerosols needing to adjust restarts from the previous forecast won't be needed, AFAIK.

For aerosol mass conservation, likely not.

For CO₂ mass conservation: depends on the application. For forward forecasts of where CO₂ gradients will be, this issue will be less relevant once someone works out how to add CO₂ observations to their runs, although there will be people who restart the meteorology from reanalysis every one to five days to avoid the expense of data assimilation. For the long-term runs used to analyze or estimate CO₂ flux fields, mass conservation is more important, so people will likely avoid DA.

As for whether those concerns are in scope for global-workflow, or should be redirected to UFS or FV3, I will leave that to the maintainers.

[jianheACM]

Daniel, thank you for checking. Looks like the rescaling approach needs further consideration.

I noticed also on slide 3 of Jian He's ppt there is notation Ptop=0. Please note FV3 models cannot have a model top at zero. Ptop should be 20 Pa. If there is a zero pressure in gfs_data.tile[1-6].nc it should have never been used in any model integration. Not sure if this has played any role in interpolation and calculation of mass.

In the ufs model source code (external_ic.F90):
if (Atm%flagstruct%external_eta) then
itoa = levp - npz + 1
Atm%ptop = ak(itoa)
Atm%ak(1:npz+1) = ak(itoa:levp+1)
Atm%bk(1:npz+1) = bk(itoa:levp+1)
call set_external_eta (Atm%ak, Atm%bk, Atm%ptop, Atm%ks)
endif
! call vertical remapping algorithms
if(is_master()) write(,) 'GFS ak(1)=', ak(1), ' ak(2)=', ak(2)
ak(1) = max(1.e-9, ak(1))

In this case, ak(1) = 1.e-9 instead of 0. The subroutine remap_scalar controls vertical interpolation with different layers of interpolation subroutines. We do not have chem mass above 20Pa in IC and the targeted model top also is 20Pa. In theory, it should not affect the total mass. But I'm not sure how exactly the interpolation functions work. There still might be something that needs attention or fix.