exports and LONG_NAMES for land water and energy balances

[gmao-rreichle]

This PR cleans up the various "spurious" (or "accounting") flux terms that are needed to reconcile what Catchment calculates and what the atmospheric turbulence expects (based on its linearization).

A matching PR for M21C (#963) was merged on 24 July 2024. The changes still need to be merged into the "develop" branch and ultimately adopted in GEOS-FP.

Additionally, another PR is needed to capture the the extensive revisions of the output (export) variable (long) names undertaken for M21C. The R21C and develop branches will keep diverging, and it will become more and more difficult to capture the considerable work done for M21C in develop.

The PR is also summarized (incl. key plots) in the following presentation:
M21C_land_budgets_20240726.pdf
M21C_land_budgets_20240726.pptx
[I apologize in advance for typos in the above slides and the discussion below. Please send corrections and suggestions to @gmao-rreichle.]

Related PR: GEOS-ESM/GEOSgcm_App#621

cc: @wmputman @rdkoster @rgelaro @sdrabenh @biljanaorescanin @elakkraoui @rlucches @lcandre2 @zhaobin74 @gmao-jkolassa

---

Discussion

The PR represents a paradigm shift in the meaning of the surface turbulent fluxes returned from the Catch GC. Previously, the fluxes were adjusted to reflect what the Turbulence GC expects (based on its linearization). However, there were two errors in the adjustment of the components of the latent heat flux (EVP[*]), resulting in inconsistent output from Catch GC. (The total latent heat flux did not match the sum of its components.) It proved impossible to adjust the EVP[*] components in a way that achieved such consistency. The new paradigm outputs the surface turbulent fluxes as calculated by the Catch GC, which makes more sense anyway. The fluxes are output along with the adjustment ("accounting") terms for the evap mass flux and the sensible heat flux (which were already available and remain unchanged). A new export variable (SPLH) for the adjustment ("accounting") term of the latent heat flux was added for completeness.

Closure of the land surface water and energy balance was successfully verified after commit 82b359f by @gmao-rreichle using a one-day AMIP-style GCM experiment:

WCHANGE = PRECTOTLAND - EVLAND - RUNSURFLAND - BASEFLOWLAND

ECHANGE = SWLAND + LWLAND - SHLAND - LHLAND - Lf*PRECSNOLAND - SPSNLAND

EVLAND = ( LHLANDINTR + LHLANDSOIL + LHLANDTRNS )/Lv + LHLANDSBLN/Ls

LHLAND = LHLANDINTR + LHLANDSOIL + LHLANDTRNS + LHLANDSBLN (corrected typo in LHLAND, 24 July 2024)

where all variables are from the "lnd" collection (with SURFACE as the source GridComp), and the physics constants Lv, Lf, and Ls are the latent heat of vaporization, fusion, and sublimation, resp. (MAPL_LATENT_HEAT_[*]). Note the absence of spurious ("accounting") terms for the surface turbulent fluxes in the land balance equations. The spurious snow energy term (SPSNLAND, formerly SPSNOW) is still needed. (Variable names are close to what we want in M21C but may not the final names. Note the renaming of the latent heat flux components from EVP[*] to LHLAND[*].)

For grid cells with land fraction close to unity (FRLAND --> 1) , it was further verified that:

EVLAND - SPEVLAND = EVAP_from_Turbulence_GC

LHLAND - SPLHLAND = LHFX, consistent with mass flux EVAP_from_Turbulence_GC

SHLAND - SPSHLAND = SH_from_Turbulence_GC

which establishes consistency between and within the "flx" and "lnd" outputs of the surface turbulent fluxes.

The PR also adds in-code documentation, including a table (in Catch GC) that lists the various argument, variable, and export names for the surface turbulent fluxes and how they are related.

For now, the additional "SPLH" export and the revisions of the LONG_NAMES of the spurious terms were only implemented in Catchment. (CatchmentCN will be addressed when changes have been accepted for Catchment.)

The PR is expected to be 0-diff for restarts. In AGCM simulations, the EVLAND, LHLAND, SHLAND, and EVP[*] diagnostics in the "lnd" collections are expected to change (paradigm shift and bug fix). In offline (land-only) simulations, the LHLAND and EVP[*] diagnostics may have roundoff differences (although early testing suggests 0-diff diagnostics).

Surface turbulent fluxes averaged over all tile types in other collections ("ocn", "gmichem", "S2S") are not expected to change. However, for the land contribution it looks like the output of the surface turbulent fluxes in these collections is an inconsistent mix of what Catchment calculates (evap, sensible) and what Turbulence GC expects (latent). This should probably be cleaned up by changing the latent heat flux variable in the respective HISTORY files.

---

The following are early research notes. These notes do not reflect the final implementation. Some of the findings and proposed changes in these notes were subsequently abandoned. See the link above for for final summary.

• RR notes on spurious and accounting terms (LINK CORRECTED ON 4 JULY 2024):

  • SP_and_ACC_src_code_20240628.pdf

• RR slides on land energy and water balances:

  • M21C_land_budgets_20240611_edited_13Jun.pdf
  • M21C_land_budgets_20240628.pdf

[gmao-rreichle]

@rdkoster, @biljanaorescanin:

I started looking into Biljana's latest 1-day AMIP simulations using the "develop" branch ("stock") and this PR's branch ("exp").
I didn't quite see what I was expecting before I ran out of time.

However, it looks like the apportioning of the EVACC term in GEOS_CatchGridComp.F90 when running in coupled land-atm (GCM) mode has two errors. The first error was that the EVACC correction was only applied for SUMEVP>0, which should already fixed in the branch. The second error seems to be a confusion of units. For example, in the following equation:

EVPICE = EVPICE - EVACC*EVPICE/SUMEV

EVPICE and SUMEV are in W/m2 and EVACC is in kg/m2/s. It looks like a units conversion with L_s (latent heat of sublimation, MAPL_ALHS) is missing. Similarly, the next three equations are missing a units conversion with L_v (latent heat of vaporization, MAPL_ALHL). In the current code, the correction is therefore essentially ignored (besides being only applied for evaporation but not for dewfall) and the EVP[*] output becomes inconsistent with the EVAP and LHLAND variables.

The error might have been prevented if the code included clear documentation of the units of the inputs and outputs of catchment(). This information is sorely needed to make the code understandable to future investigators.

If the above makes any sense, please modify the branch accordingly and re-run the 1-day AMIP experiments.

@biljanaorescanin: To be safe, it would be best to switch to a restart date in mid-February. We tend to have the biggest problems with the snow terms, and the snow-covered land area is larger in Feb than in Apr, so any issues with snow will be easier to detect.

[gmao-rreichle]

Thanks, @biljanaorescanin, for making the proposed changes and running another day in Feb.

For this PR's branch at d95a4a1, the following energy balance:

ECHANGE = SWLAND + LWLAND - SHLAND 
          - LHLAND 
          - Lf*PRECSNO - SPLAND - SPSNOW - SPLH

closes within +/- 0.2 W/m2 for all grid cells, with a std-dev across all grid cells of ~0.02 W/m2.
Moreover, we now have

EVLAND = (EVPINTR + EVPSOIL + EVPTRNS)/Lv + EVPSBLN/Ls

to within 1e-3 mm/d. [CORRECTED ON 12 JULY 2024, ORIGINAL NUMBER OF 1e-5 mm/d WAS A TYPO.]
Alas, we still only have

LHLAND = EVPINTR + EVPSBLN + EVPSOIL + EVPTRNS

to within +/- 2 W/m2. The mismatch seems to occur only when FRSNOW>0, but I'm not sure. This requires more investigation.

For reference, in "develop" the following energy balance:

ECHANGE = SWLAND + LWLAND - SHLAND 
          - EVPINTR - EVPSBLN - EVPSOIL - EVPTRNS
          - Lf*PRECSNO - SPLAND - SPSNOW 

is equivalent to the branch's balance above (i.e., closes to within +/- 0.2 W/m2 for all grid cells, with a std-dev across all grid cells of ~0.02 W/m2).

However, for "develop":

EVLAND = (EVPINTR + EVPSOIL + EVPTRNS)/Lv + EVPSBLN/Ls

does not match (deviates up to 10 mm/d) and

LHLAND = EVPINTR + EVPSBLN + EVPSOIL + EVPTRNS

does not match (deviates up to 400 W/m2).

cc: @rdkoster

[gmao-rreichle]

0-diff tests for GEOSldas successfully completed by @biljanaorescanin on 8 July 2024

[biljanaorescanin]

@gmao-rreichle last commit is @rdkoster suggestion how to change units. I will send run output once done in email.

[gmao-rreichle]

@gmao-rreichle last commit is @rdkoster suggestion how to change units. I will send run output once done in email.

@rdkoster, @biljanaorescanin: Thanks again for implementing the change. At first glance it does make more sense to do the corrections in this way. (As an aside, in the new corrections, we should be able to just use EVLAND instead of computing SUMEV; likewise, we could probably have used LHOUT instead of computing SUMEV in the old corrections).

Alas, the change does not solve the remaining issues. To recap, the issues are:

1. LHLAND=EVPSUM is true only to within +/- 2 W/m2 at a subset of tiles.

2. When using LHLAND in the energy balance, the balance only closes to within +/- 0.2 W/m2 at a subset of tiles. (The balance is worse -- to within +/- 2 W/m2 -- if we use EVPSUM instead of LHLAND.)

The latest commit only changes EVPSUM and not LHLAND, the energy balance using LHLAND is unchanged. When using the revised EVPSUM in the energy balance, the balance is no better than with old EVPSUM (to within +/-2 W/m2). Also, with the revised EVPSUM, the equality EVPSUM=LHLAND is still true only to within +/- 2 W/m2.

It's probably still true that the lack of closure only happens where FRSNOW>0, although there are tiles with FRSNOW>0 that do have closure. But I ran out of time to look into this further. I might also have messed up my plots in the rush to take a quick look. We'll have to revisit next week.

[gmao-rreichle]

@zhaobin74 : I changed three long names of CICE longwave exports in commit feb6d1d
These changes seem to have been missed in #764, but I'm not sure. Can you please double-check and let me know?
Note that the present PR is on top of develop. I have another PR in the works that is on top of R21C. The latter is urgent. I would greatly appreciate it if you could take a look soon. Thanks!

[zhaobin74]

@zhaobin74 : I changed three long names of CICE longwave exports in commit feb6d1d These changes seem to have been missed in #764, but I'm not sure. Can you please double-check and let me know? Note that the present PR is on top of develop. I have another PR in the works that is on top of R21C. The latter is urgent. I would greatly appreciate it if you could take a look soon. Thanks!

@gmao-rreichle, just back from a long leave. I missed changes to these variables in #764. Thanks for adding them.