Bug in heat capacity equation and soil density and soil thermal conductivity

[penguian]

keyword_maygit resolution_fixed type_defect | by m.decker@unsw.edu.au

---

CABLE has two bugs related to soil thermal properties, found by Anne during SOIL-MIP

BUG #1: Uses bulk density when actual density is the correct parameter.

BUG #2:The min/max/if/else statements create nonphysical decreasing conductivity as moisture increases, contrasting what I believe is the actual equation

FIXED AS PART 0F TICKET #191 BUG #1

EQN 1: Csoil,i = Δ zi [ρmrl (1 - Θsat )Cvmrl + ρmrl Θliq Cvliq + ρmrl Θice *Cvice

ρmrl is the actual density of the minerals (universally taken as 2700.0 kg/m3 (though really +-100.0)). ρliq and #961ice are the actual density of liquid and ice water

CABLE mistakenly reads the bulk density of the soil and uses this value in (1). The density of the mineral material itself is the correct parameter. def_soil_params.txt lists values of soil density (actually ρbulk not ρ{mrl}) around 1300.0-1700.0.

These are unphysically small for ρmrl but in the ball park of observed ρbulk.

The relationship bewteen actual and bulk soil density is:

EQN 2: ρ{bulk} = ρ{mrl} * (1 - Θsat)

So it seems CABLE reads in the bulk density, then incorrectly multiplies the bulk density by the total volume (1 - Θsat.

However Eqn 1 as coded in CABLE introduces a max statement into Eqn 1:

EQN 3: Csoil,i = MAX [ Δ zi [ρmrl (1 - Θsat )Cvmrl + ρmrl Θliq Cvliq + ρmrl Θice Cvice , ρmrl C{v}mrl ]

giving the correct values for C{soil,i} when total moisture goes to zero, but neglecting the initial increase in Csoil,i when wetting occurs.

BUG #2

Well the equation to calculate thermal conductivity as a function of soil liquid and ice has no relation to the 2006 documentation, is not from a published method, and was made to match some curves.

Actual equation is very complex (at end of message), with numerous min/max/if conditions, but removing these and simplifying the log raised to exp nonsense gives:

EQN 3: Λsoil = Λdry,soil 1.41 60.0Θliq 250Θice Θliq-1

Note: I made the simplifications and this is not in the code. Figures below are made with actual code. The above simplified code produces physically reasonable conductivity values (increases with liq or ice)

Figures attached.

---

Issue migrated from trac:176 at 2023-11-27 11:23:56 +1100

[penguian]

m.decker@unsw.edu.au _uploaded file summary_two_bugs.pdf (297.8 KiB)_

Full description including figures

[penguian]

m.decker@unsw.edu.au changed owner from bep599 to somebody

[penguian]

m.decker@unsw.edu.au changed component from MPI to biogeophys

[penguian]

m.decker@unsw.edu.au commented

---

Ticket #191 includes the code changes needed to fix the problematic heat capacity and conductivity.

A new cable_user switch has been added: cable_user%soil_thermal_fix. When false, the model behaviour is unchanged, but when true the heat capacity and conductivity are correctly calculated.

Heat capacity fix The heat capacity calculation (in soilsnow) is changed to

gammmzz = max(soil%heat_cap_lower_limit,
              ( 1.0 - soil%ssat ) * soil%css * soil%rhosoil * &
              + soil%ssat * ( wblfsp * cswat * C%density_liq + &
              ssnow%wbfice(:,k) * csice * C%density_liq * 0.9 ) )*soil%zse

where soil%heat_cap_lower_limit = 0.01 if cable_user%soil_thermal_fix=true, and

soil%heat_cap_lower_limit = soil%css*soil%rhosoil

when it is false.

Additionally, the soil density (soil%rhosoil) is set to 2700.0 (the density value of the soil grains, NOT the bulk density).

Conductivity fix

The default code for calculating the thermal conductivity has been moved into a function called old_soil_conductivity(ssnow,soil). Although in a new function the code is exactly the same as before. Setting cable_user%soil_thermal_fix causes the conductivity to be found using total_soil_conductivity(ssnow,soil) (currently called in stempv and GWstempv) or can be changed (in code) to call sli_soil_cond(ssnow,soil).

total_soil_conductivity(ssnow,soil) uses the formulas from Peters-Lidard,. 1998). It does not require any new parameters. The dry soil conductivity is found from the bulk soil density.

sli_soil_cond(ssnow,soil) uses the formulas from SLI (Haverd and Cuntz 2010, Cuntz and Haverd 2017). It has parameters that have been added as soil%Kthermconst{one,two,three} that are calculated during initialization.

The choice between the above is not on a switch but could be. The final code should use total_soil_conductivity.

[penguian]

@jxs599@nci.org.au changed status from new to closed

[penguian]

@jxs599@nci.org.au changed resolution from ` tofixed`

[penguian]

@jxs599@nci.org.au edited the issue description

[penguian]

@jxs599@nci.org.au changed milestone from ` to1. Closed`

[penguian]

@ccc561@nci.org.au set keywords to maygit