Hillslope spinup and calibration

[wwieder]

In #147 I started a spinup using @swensosc's hillslope branch to start evaluating carbon cycle difference with the hillslope model turned on. Expected changes were higher soil C pools where higher water table causes lower values for the O-scalar, which would increase soil C pools (this seems accurate, although not shown).

It seems like our latest 'flat' f09 case is #131, which doesn't have all the latest PPE calibrations, but hopefully it isn't far off. The hillslope case also is using the two most dominant PFTs and no-crop while the flat case is full PFT crop simulation. Results are all for 1850 land use, CO2, etc.

After 180 years of AD spinup we're seeing really big differences in savannah LAIs (here #147-131), left panel. The 'raw data' from the baseline case (131) is on the right

What could be causing this? A few poorly formed thoughts, which are likely interacting:

• High water table or hillslope convergence regions increases LAI.
  • Column level output would be interesting to see if lowland columns are especially productive.
• PFT calibrations?
  • Most of the region is PFTs 6 & 14, with C4 grasses being especially dominant in Africa.
  • C4 grasses, in the base case we struggled some with C4 grass survival are they just surviving better with although regions with the largest LAI increases already have high LAI?
  • Tropical deciduous trees are really going nuts with hillslope hydrology for some reason (or make up a larger fraction of the gridcell with the 2 PFT setup)?
  • It seems like it would be helpful to diagnose which PFTs are most responsible for the change?
• Reduced fire, again potentially from wetter fuels with the hillslope model?

I'm going to get some of the column and PFT level output to start diagnosing this
@swensosc will work on creating a regional sandbox to start investigating

@linniahawkins, @dlawrenncar, @lifang0209, @slevis-lmwg, @olyson curious to hear your thoughts here.

[wwieder]

First pass looking at a hillslope vs. 'flat' case are here. A few notes, the hillslope case is a BGC-noCrop case that is not spun up, this is a 10 year simulation from year 180 of the AD spinup. The flat case is from the first decades of pSASU in #131, which is BGC-crop using an older parameter file which should have lower productivity in tropical and boreal forests.

LAI and GPP are way up in the hillslope case (global GPP is 130 vs 100 PgC /y in hillsope and flat cases, respectively).
Notably burned area is much reduced in the hillslope case (1 vs. 9 millon km2), likely reflecting wetter conditions in the hillslope results.

Some of the biggest differences are in savannah regions. Here's a look at annual means for the two cases in southern Africa.

[wwieder]

Differences in burned area are pretty stark. Maybe we need to recalibrate the fire model before going too much further?

[wwieder]

Quick update on a few thing.

I ran a test to see if running with full PFTs would influence much. It does not and reduces throughput

Below is ELAI in year 8 from a

• full PFT (42 ypd), left;
• ndomPFT=2 (65 ypd, Sean's default configuration), center; and
• non-hillslope, full PFT case (2 degree, but using the same model tag and parameter file)

All color bars are identical

Hillslope runs look to have higher LAI in savanahs and lower LAIs in tropical forests, compared to the flat configuration.
Adding full PFTs may reduce some high LAI biases (e.g. SE US), relative to the 2 PFT case but not enough to warrant the slower throughput.

[wwieder]

@lifang0209, we're trying a CRU-JRA forced run using the hillslope hydrology model and are seeing very high LAI and low burned area in Savanah regions. I'm assuming this is because of higher water table and wetter surface soils in the hillslope runs (see post above)?

I'm wondering if you can suggest any change to the file mode we can explore here?

The sample history output for the hillslope case is here:
/glade/derecho/scratch/wwieder/archive/ctsm5.4.022_BGC_hillslope_f09_147_test1HH

A corresponding 'flat' case is here
/glade/derecho/scratch/wwieder/archive/ctsm5.4.022_BGC_hillslope_f09_147_test1flat

Neither of these cases are spun up, they both are from ~180 y of AD spinup.
We're also using a 16 pft BGC no crop setup with only 2 dominant PFTs in each grid cell.

Let me know if you have ideas / suggestions here, or if finding a time to talk would be helpful. Thank you.

[lifang0209]

@wwieder Are these offline simulations using the fire parameters tuned for the CESM coupled case? And is the hillslope in the CLM6 default?
After turning on hillslope, the fuels become wetter and burned area decreases. However, it’s really surprising that activating the hillslope module would drive burned area in tropical savannas down to nearly zero.
I’ll get back to you tomorrow or the day after. There’s an important deadline tomorrow.

[slevis-lmwg]

@lifang0209 I will answer a couple of question for which I know the answers:
Yes, these are offline simulations. Fire parameters are NOT for the CESM coupled model. These are the paramfiles used, the first one for the hillslope case and the second one for the "flat" case:

<  paramfile = '/glade/campaign/cesm/cesmdata/inputdata/lnd/clm2/paramdata/ctsm60_params.c260204.no_nan_fill.nc'
>  paramfile = '/glade/campaign/cesm/cesmdata/inputdata/lnd/clm2/paramdata/ctsm60_params.c251211.nc'

You can see the lnd_in files here:

< /glade/derecho/scratch/wwieder/ctsm5.4.022_BGC_hillslope_f09_147_test1HH/run
> /glade/derecho/scratch/wwieder/ctsm5.4.004_bgc_f09_131_pSASU/run

[lifang0209]

(1) After incorporating Hillslope, the simulated burned-area fraction (%; Fig. 1) becomes unrealistically low. The global total burned area drops from 837 Mha/yr to 93 Mha/yr.

(2) This underestimation in Hillslope case is driven by much higher root-zone soil wetness (BTRAN2; 0–1), which is a factor in fuel wetness. The increase is substantial (not marginal). Fig. 2 shows monthly BTRAN2 during the observed peak fire month.

(3) Have you evaluated the soil wetness/soil moisture in the Hillslope case? If the higher soil wetness is reasonable, I can tune rswf_max and rswf_min to match the new soil wetness range.

(4) Alternatively, the Hillslope developer could check whether the physical meaning of any variables used in calculating BTRAN2 (see below) changed after Hillslope was introduced.

do j = 1,nlevgrnd
do fp = 1, num_exposedvegp
p = filter_exposedvegp(fp)
c = patch%column(p)
l = patch%landunit(p)
s_node = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8)
call soil_water_retention_curve%soil_suction(c, j, s_node, soilstate_inst, smp_node_lf)
smp_node_lf = max(smpsc(patch%itype(p)), smp_node_lf)
btran2(p) = btran2(p) +rootfr(p,j)*max(0._r8,min((smp_node_lf - smpsc(patch%itype(p))) / &
(smpso(patch%itype(p)) - smpsc(patch%itype(p))), 1._r8))
end do
end do

[wwieder]

Thanks for this summary @lifang0209. I wonder how hard it would be for you to tune rswf_max and rswf_min (option 3 above) given these new conditions? Given the time constraints we're under this seems like the most expedient option, but we'll have lots of exciting science to consider down the road!

[lifang0209]

@wwieder (1) I can send you the retuned rswf_max and rswf_min next Monday or Tuesday. I need to run a few experiments to determine the updated values.
Could you please confirm that I should create the clone like this?
cd /glade/u/home/wwieder/wetland_update/cime/scripts
./create_clone --clone /glade/work/wwieder/wiemip/ctsm5.4.022_BGC_hillslope_f09_147_test1HH /glade/work/fangli/hillslope

(2) I’m curious why root-zone soil moisture increases so much when hillslope processes are introduced, and whether this increase is physically realistic. In China’s ESMs (BCC and CAS), lower rswf_max and rswf_min than in CLM5 imply a drier root zone, yet adding the hillslope parameterization to CLM makes it much wetter.

[wwieder]

Sorry, for my delay yes to question # 1, that's how you'd create the clone, just add the --case flag

It looks like I've updated my wetland_update branch to ctsm5.4.024, but hopefully this works as intended for you. Let me know if you run into issues with building your hillslope case.

I'd assume the hillslope parameterization will result in wetter soils in lowland columns, but that upland columns in the hillslope should be similar to single column represented 'flat' case. There's a lot of interesting science we'll get to start doing with the hillslope configuration, but for now we just need something that looks plausible for LAI and burned area.

[wwieder]

@lifang0209 I wonder if you've had a chance to look into suggested parameter values for rswf_max and rswf_min with the hillslope model yet?

[lifang0209]

I’m still checking the root-zone soil moisture. The root-zone soil wetness (btran2) looks too high in forest areas in hillslope case. I’ll get back to you soon on whether adjusting rswf_max and rswf_min is sufficient.

[swensosc]

It is important to note that the hillslope model has multiple columns per gridcell. So when gridcell values are plotted, they are showing an average across the columns. In many places that are colored red in the figure above, wetlands are present. So these gridcells will likely have a very wet lowland column, but drier upland columns. So the gridcell value may be weighted towards the wetter lowland column. Here is a recent comparison of a hillslope simulation and observations of wetlands:

[wwieder]

I’m still checking the root-zone soil moisture. The root-zone soil wetness (btran2) looks too high in forest areas in hillslope case. I’ll get back to you soon on whether adjusting rswf_max and rswf_min is sufficient.

Thanks Fang. Just to clarify we're expecting the hillslope hydrology configuration to have different soil wetness and require different fire parameters. This is an active research area, and the parameters you suggest won't be used in 'off the shelf' CLM6 / CESM3 simulations, they will be specific for configurations of the model using hillslope hydrology.

[lifang0209]

@swensosc Thanks for the explanation. That's very helpful!

I think the issue on fire simulations may be related to how the averaged root-zone soil wetness is correctly calculated for a PFT in a vegetated column when the hillslope model is active and lowland/upland columns are introduced.
Could you take a look at the code below, which calculates btran2 across different PFTs in a vegetated column, and see whether any modifications are needed when hillslope is active?

do fp = 1, num_exposedvegp
p = filter_exposedvegp(fp)
btran2(p) = 0._r8
end do

do j = 1,nlevgrnd
do fp = 1, num_exposedvegp
p = filter_exposedvegp(fp)
c = patch%column(p)
s_node = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8)
btran2(p) = btran2(p) + rootfr(p,j)*s_node
end do
end do

do fp = 1, num_exposedvegp
p = filter_exposedvegp(fp)
if (btran2(p) > 1._r8) then
btran2(p) = 1._r8
end if
end do

[wwieder]

This is an interesting question, @lifang0209, because grid cell averaged differences in BTRANMN that's used by PHS have a different spatial pattern (notably drier in some savannah regions) than what we're seeing in differences for BTRAN2 (wetter everywhere)

more diagnostics here

[swensosc]

I think the btran2 calculation is still correct, because it occurs separately for each pft on each column. However, when they are averaged together, it will not show those differences between columns. This figure shows soil water in layer 1 for four columns (0->lowland, 3->upland). If you compare 0 to 3, you can see that the upland is often 1/2 the value of the lowland. I expect in the hillslope model, uplands in a gridcell may burn but the lowland may not, which would affect the burned area. As a simple example, imagine if we had two columns and wanted to get a burned area of 0.5. If the lowland was was fairly wet, it may not burn at all, so the upland would need to burn completely to have a gridcell burned area of 0.5.

[wwieder]

Thanks for looking into this, Sean. @lifang0209 do you think you'll be able to provide suggested parameters that generate better (but maybe not perfect) burned area estimates pretty quickly? This is something we'll need pretty quickly if we're going to try using the hillslope configuration of CLM for WIE-MIP.

[lifang0209]

@wwieder I'll work on this ASAP. However, I suspect the issue may not be just about parameters. Some column-level fire-related variables may not be correctly calculated when hillslope is active. It would be very helpful if I could get some assistance from the hillslope developer to look into this together.

[lifang0209]

@swensosc Thanks, Sean!
In CLM, burned area fraction (farea_burned) is calculated at the column level. When calculating farea_burned for a natural vegetated (including bare soil) column, we first need to compute the column-averaged btran_col based on PFT-level btran2 as:

do fc = 1,num_soilc
c = filter_soilc(fc)
btran_col(c) = 0._r8
end do
do fp = 1, num_exposedvegp
p = filter_exposedvegp(fp)
c = patch%column(p)
! For non-crop -- natural vegetation and bare-soil
if( patch%itype(p) < nc3crop .and. cropf_col(c) < 1.0_r8 )then
btran_col(c) = btran_col(c)+max(0._r8, min(1._r8, (btran2(p)-rswf_min(patch%itype(p))) &
/(rswf_max(patch%itype(p))-rswf_min(patch%itype(p)))))*patch%wtcol(p)
end if
end do

Sean, could you take a look at whether any modifications are needed when the hillslope model is active? It seems that we are only using the lowland column level (which is very wet).

I suspect the column indexing is disrupted when the hillslope model is active. I set peat fires to zero in grid cells where the tropical forest area fraction exceeds 0.8 ((trotr1_col(c) + trotr2_col(c)) * col%wtgcell(c) > 0.8_r8). However, when hillslope is active, peat fires are unexpectedly greater than zero in these grid cells.

I have a related question: when hillslope is active, for a grid cell that originally has 2 columns (e.g., a natural vegetated column and a rice column), does the number of columns increase to 4 (i.e., lowland and upland columns for the vegetated and rice columns, respectively)?

[swensosc]

hi, fang,
Currently, hillslope only operates on the vegetated landunit. So crop column number does not change, just the nat veg columns. I think the btran_col calculation looks okay because it is identifying the correct c for each p. But the forest fraction will likely be too small now because col%wtgcell will be much smaller with multiple columns. With hillslope, you would need to either 1) sum over columns (i.e. sum(trotr1+tottr2)*wt > 0.8) or 2) calculate the peat fire activation separately for each column.

[swensosc]

thinking about this a little more, why is col%wtgcell(c) needed? If you are calculating a column level variable, i.e. baf_peatf(c), is it necessary to know how much of a gridcell that column occupies? For example, if a gridcell was 0.5 crop and 0.5 BET (trotr1_col(c) = 1) doesn't that indicate a closed canopy forest on that column?

[lifang0209]

@swensosc Thanks. The reason I multiply by col%wtgcell(c) is that observations suggest there are no peat fires when tropical tree fraction exceeds 80% at the grid-cell level, so I need to convert the column-level fraction to a grid-level area fraction.

I think I may have found the reason why burned area becomes so small when hillslope is active. Previously, each grid cell had only one natural vegetated column, so multiplying by col%wtgcell(c) correctly converted column-level fraction to grid-level area fraction and it's simple. With hillslope on, that natural vegetated column is split into upland and lowland columns, so col%wtgcell(c) for each one is smaller. The same issue likely affects lfwt(c) in the non-peat burned area calculation outside tropical closed forests and croplands.

I’ll submit two simulations to test this.

[swensosc]

hi, fang,
have you had time to look at this issue? I did a short simulation in which I removed the col%wtgcell from the troptr_ if-statement, and also changed col%wtgcell to col%wtcol in the lfwt calculation. The farea_burned magnitudes then looked similar to those in a non-hillslope simulation. This plot shows my hillslope run with these changes. Each panel is a different column (the columns all look fairly similar). I multiplied the history output by (30x24x3600 seconds per month).

[lifang0209]

@wwieder Hi Will, could you re-run the spin-up using the corrected F90 files in /glade/u/home/fangli/will2026/hh-test/SourceMods/src.clm ?

I added tgswt to represent the gridcell fractional coverage of tree, grass, and shrub PFTs. This is meant to capture the continuity of natural vegetation within a grid cell: if the fractions of urban, cropland, bare soil, or lake are larger, then tgswt is smaller, and fuel continuity is reduced.

When hillslope is off, there is only one natural vegetation column in a grid cell, so lfwt(c) (sum of tree, grass, shrub's patch%wtgcell(p)) correctly represents tgswt(g) without introducing a new variable name. However, when hillslope is on, that is no longer correct.

After you get a 1850 spin-up, we can see whether we need to tune the thresholds for btran2.

@swensosc Thanks, Sean! but I need the gridcell fraction.
In addition, I’m curious how many columns are in the natural vegetation landunit now. Is it always 4, or does it depend on the topography in a grid cell? In the figure you showed above, if col0 is lowland and col3 is upland, what are col1 and col2?

[swensosc]

hi, fang,
the number of columns is not always the same. It depends on the input dataset, and it can be different for each gridcell. The plot I showed was from a simulation where there are actually 16 columns. Here is figure showing how columns 1-4 are connected:

[lifang0209]

Hi Sean, thank! I understand now.

[swensosc]

regarding the wtgcell: I was just removing it to test whether that was the issue. I am glad you have created a fix that will work for both hillslope and non-hillslope. I am testing now in a hillslope run... Here is the output of a 5 year run using fang's mods:

[wwieder]

I'm also testing this in a non-hillslope run in #161