Add extended CLM snow parameterization options Feature/clm snow (parflow#698)

## Summary
Extends CLM snow parameterizations with new rain-snow partitioning
methods, SZA-based melt damping, and configurable fractional snow cover.
**New Features:**
- **Dai (2008) rain-snow partitioning**: Hyperbolic tangent function
fitted to 30 years of global observations
- **Jennings (2018) rain-snow partitioning**: Bivariate logistic
regression using temperature and relative humidity
- **SZA-based snow melt damping**: Reduces melt energy at high solar
zenith angles where CLM underestimates albedo
- **Configurable fractional snow cover**: Adjustable roughness parameter
for frac_sno calculation
- **Configurable temperature thresholds**: SnowTCrit, SnowTLow,
SnowTHigh, SnowTransitionWidth for existing methods
**New Keys (16 total):**
- `Solver.CLM.SnowTCrit`, `SnowTLow`, `SnowTHigh`, `SnowTransitionWidth`
- `Solver.CLM.DaiCoeffA`, `DaiCoeffB`, `DaiCoeffC`, `DaiCoeffD`
- `Solver.CLM.JenningsCoeffA`, `JenningsCoeffB`, `JenningsCoeffG`
- `Solver.CLM.SZASnowDamping`, `SZADampingCoszenRef`,
`SZADampingCoszenMin`
- `Solver.CLM.FracSnoScheme`, `FracSnoRoughness`
**References:**
- Dai (2008) GRL doi:10.1029/2008GL033295
- Jennings et al. (2018) Nature Communications
doi:10.1038/s41467-018-03629-7
- Dang et al. (2019) The Cryosphere doi:10.5194/tc-13-2325-2019
## Test plan
- [x] `clm_snow_dai` - Dai sigmoidal partitioning method
- [x] `clm_snow_jennings` - Jennings bivariate logistic method
- [x] `clm_snow_sza_damping` - SZA-based melt reduction
- [x] All existing CLM snow tests pass
Run tests with: `ctest -R clm_snow`
  🤖 Generated with [Claude Code](https://claude.ai/code)

---------

Co-authored-by: Reed <rmaxwell@mines.edu>
Co-authored-by: Claude <noreply@anthropic.com>

Author: 3 people

docs/user_manual/keys.rst

@@ -6357,23 +6357,77 @@ Reference: :cite:p:`Wang2019`.
 
 *string* **Solver.CLM.SnowPartition** CLM Selects the method for
 partitioning precipitation into rain and snow. The valid types for
-this key are **CLM**, **WetbulbThreshold**, **WetbulbLinear**.
+this key are **CLM**, **WetbulbThreshold**, **WetbulbLinear**, **Dai**, **Jennings**.
 
 **CLM**:
    Standard air temperature threshold with linear transition (default).
+   Uses configurable SnowTLow and SnowTHigh thresholds.
 
 **WetbulbThreshold**:
    Sharp threshold at wet-bulb temperature. Better for dry mountain climates.
 
 **WetbulbLinear**:
    Linear transition around wet-bulb temperature threshold.
+   Uses configurable SnowTransitionWidth for transition zone.
+
+**Dai**:
+   Sigmoidal function of air temperature from :cite:p:`Dai2008`.
+   Uses configurable DaiCoeffA/B/C/D coefficients.
+
+**Jennings**:
+   Bivariate logistic regression with air temperature and relative humidity
+   from :cite:p:`Jennings2018`. Uses configurable JenningsCoeffA/B/G coefficients.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.SnowPartition "Dai"         ## TCL syntax
+      <runname>.Solver.CLM.SnowPartition = "Dai"   ## Python syntax
+
+*double* **Solver.CLM.SnowTCrit** 2.5 Initial classification threshold
+above freezing (K) for determining precipitation type in drv_getforce.
+If air temperature exceeds tfrz + SnowTCrit, precipitation is initially
+classified as rain. Default 2.5 K matches the hardcoded CLM value.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.SnowTCrit 2.5         ## TCL syntax
+      <runname>.Solver.CLM.SnowTCrit = 2.5   ## Python syntax
+
+*double* **Solver.CLM.SnowTLow** 273.16 CLM method lower temperature
+threshold (K) below which all precipitation is snow. Default 273.16 K (freezing).
 
 .. container:: list
 
    ::
 
-      pfset Solver.CLM.SnowPartition "WetbulbThreshold"         ## TCL syntax
-      <runname>.Solver.CLM.SnowPartition = "WetbulbThreshold"   ## Python syntax
+      pfset Solver.CLM.SnowTLow 273.16         ## TCL syntax
+      <runname>.Solver.CLM.SnowTLow = 273.16   ## Python syntax
+
+*double* **Solver.CLM.SnowTHigh** 275.16 CLM method upper temperature
+threshold (K) above which the liquid fraction reaches maximum (40%).
+Default 275.16 K (tfrz + 2).
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.SnowTHigh 275.16         ## TCL syntax
+      <runname>.Solver.CLM.SnowTHigh = 275.16   ## Python syntax
+
+*double* **Solver.CLM.SnowTransitionWidth** 1.0 WetbulbLinear method
+half-width (K) of transition zone. Transition spans threshold +/- this value.
+Default 1.0 K for 2K total range.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.SnowTransitionWidth 1.0         ## TCL syntax
+      <runname>.Solver.CLM.SnowTransitionWidth = 1.0   ## Python syntax
 
 *double* **Solver.CLM.WetbulbThreshold** 274.15 Threshold temperature in
 Kelvin for wetbulb partitioning methods. Default 274.15 K (1°C). Only
@@ -6387,6 +6441,48 @@ used when ``Solver.CLM.SnowPartition`` is ``WetbulbThreshold`` or
       pfset Solver.CLM.WetbulbThreshold 274.15         ## TCL syntax
       <runname>.Solver.CLM.WetbulbThreshold = 274.15   ## Python syntax
 
+**Dai Coefficients**
+
+The Dai (2008) method uses a hyperbolic tangent function fitted to 30 years
+of global weather station observations. Reference: :cite:p:`Dai2008`.
+
+Formula: F(%) = a * [tanh(b*(T-c)) - d], where T is in Celsius.
+Defaults are from Table 1a (Land, annual) of Dai (2008).
+
+*double* **Solver.CLM.DaiCoeffA** -48.2292 Dai coefficient a (scaling factor, negative).
+
+*double* **Solver.CLM.DaiCoeffB** 0.7205 Dai coefficient b (slope at half-frequency T).
+
+*double* **Solver.CLM.DaiCoeffC** 1.1662 Dai coefficient c (half-frequency temperature in C,
+where snow probability is ~50%).
+
+*double* **Solver.CLM.DaiCoeffD** 1.0223 Dai coefficient d (asymmetry parameter, ~1.0).
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.DaiCoeffA -48.2292         ## TCL syntax
+      <runname>.Solver.CLM.DaiCoeffA = -48.2292   ## Python syntax
+
+**Jennings Coefficients**
+
+The Jennings (2018) method uses bivariate logistic regression:
+psnow = 1 / (1 + exp(a + b*T + g*RH)), where T is in Celsius and RH in percent.
+
+*double* **Solver.CLM.JenningsCoeffA** -10.04 Jennings intercept coefficient.
+
+*double* **Solver.CLM.JenningsCoeffB** 1.41 Jennings temperature coefficient.
+
+*double* **Solver.CLM.JenningsCoeffG** 0.09 Jennings relative humidity coefficient.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.JenningsCoeffA -10.04         ## TCL syntax
+      <runname>.Solver.CLM.JenningsCoeffA = -10.04   ## Python syntax
+
 **Thin Snow Damping**
 
 Thin early-season snowpacks can experience spurious melt due to warm
@@ -6415,6 +6511,52 @@ threshold in mm below which thin snow damping applies.
       pfset Solver.CLM.ThinSnowThreshold 50.0         ## TCL syntax
       <runname>.Solver.CLM.ThinSnowThreshold = 50.0   ## Python syntax
 
+**SZA-Based Snow Damping**
+
+CLM's narrowband snow optical parameters are computed assuming a solar
+zenith angle (SZA) of 60 degrees. At higher zenith angles, actual snow
+albedo is higher than CLM assumes, meaning less energy should be available
+for melt. This is particularly relevant for high-latitude sites during
+accumulation season. Reference: :cite:p:`Dang2019`.
+
+*double* **Solver.CLM.SZASnowDamping** 1.0 Fraction of melt energy
+retained at high solar zenith angles. A value of 1.0 means no damping
+(default, disabled). A value of 0.8 means 20% energy reduction at high SZA.
+Damping varies linearly with cosine of zenith angle between the reference
+and minimum thresholds.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.SZASnowDamping 0.8         ## TCL syntax
+      <runname>.Solver.CLM.SZASnowDamping = 0.8   ## Python syntax
+
+*double* **Solver.CLM.SZADampingCoszenRef** 0.5 Reference cosine of solar
+zenith angle below which SZA damping applies. Default 0.5 corresponds to
+SZA of 60 degrees (matching CLM's assumption for optical parameters).
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.SZADampingCoszenRef 0.5         ## TCL syntax
+      <runname>.Solver.CLM.SZADampingCoszenRef = 0.5   ## Python syntax
+
+*double* **Solver.CLM.SZADampingCoszenMin** 0.1 Cosine of solar zenith
+angle at which maximum SZA damping applies. Default 0.1 corresponds to
+SZA of approximately 84 degrees. Must be less than SZADampingCoszenRef.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.SZADampingCoszenMin 0.1         ## TCL syntax
+      <runname>.Solver.CLM.SZADampingCoszenMin = 0.1   ## Python syntax
+
+Note: Both thin snow damping and SZA damping can be enabled simultaneously.
+When both are active, they combine multiplicatively.
+
 **Albedo Schemes**
 
 Snow albedo controls the radiation balance and strongly influences melt
@@ -6512,6 +6654,38 @@ snow faster.
       pfset Solver.CLM.AlbedoThawA 0.82         ## TCL syntax
       <runname>.Solver.CLM.AlbedoThawA = 0.82   ## Python syntax
 
+**Fractional Snow Cover**
+
+The fractional snow covered area (frac_sno) affects surface energy balance
+by weighting snow and bare ground contributions. The CLM formulation uses
+a tanh-like relationship with snow depth and surface roughness.
+
+*string* **Solver.CLM.FracSnoScheme** CLM Selects the fractional snow cover
+calculation method. Currently only CLM is available; extensible for future
+formulations.
+
+**CLM**:
+   Standard formulation: frac_sno = snowdp / (10*roughness + snowdp)
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.FracSnoScheme "CLM"         ## TCL syntax
+      <runname>.Solver.CLM.FracSnoScheme = "CLM"   ## Python syntax
+
+*double* **Solver.CLM.FracSnoRoughness** 0.01 Roughness length scale for
+fractional snow cover calculation [m]. Default 0.01 m matches CLM's zlnd
+parameter for backward compatibility. Larger values reduce snow cover
+fraction for a given snow depth.
+
+.. container:: list
+
+   ::
+
+      pfset Solver.CLM.FracSnoRoughness 0.01         ## TCL syntax
+      <runname>.Solver.CLM.FracSnoRoughness = 0.01   ## Python syntax
+
 
 .. _ParFlow NetCDF4 Parallel I/O:
 

docs/user_manual/refs.bib

@@ -256,6 +256,16 @@ @article{Cui14
 doi = {10.1016/j.cageo.2014.05.005},
 }
 
+@article{Dai2008,
+author = {Dai, A.},
+journal = {Geophysical Research Letters},
+pages = {L12802},
+title = {Temperature and pressure dependence of the rain-snow phase transition over land and ocean},
+volume = {35},
+year = {2008},
+doi = {10.1029/2008GL033295},
+}
+
 @article{Dai03,
 author = {Dai, Y. and X. Zeng and R. E. Dickinson and I. Baker and G. B. Bonan and M. G. Bosilovich and A. S. Denning and P. A. Dirmeyer and P. R., G. Niu and K. W. Oleson and C. A. Schlosser and Z. L. Yang},
 journal = {The Bulletin of the American Meteorological Society},
@@ -277,6 +287,16 @@ @article{Danesh-Yazdi2018
 url = {https://doi.org/10.1002/hyp.11481},
 }
 
+@article{Dang2019,
+author = {Dang, C. and Zender, C. S. and Flanner, M. G.},
+journal = {The Cryosphere},
+pages = {2325--2343},
+title = {Intercomparison and improvement of two-stream shortwave radiative transfer schemes in Earth system models for a unified treatment of cryospheric surfaces},
+volume = {13},
+year = {2019},
+doi = {10.5194/tc-13-2325-2019},
+}
+
 @article{DMC10,
 author = {Daniels, M. H. and Maxwell, R. M. and Chow, F. K.},
 journal = {Journal of Hydrologic Engineering},
@@ -557,6 +577,16 @@ @article{Jefferson2017
 url = {https://doi.org/10.1175/jhm-d-16-0053.1},
 }
 
+@article{Jennings2018,
+author = {Jennings, K. S. and Winchell, T. S. and Livneh, B. and Molotch, N. P.},
+journal = {Nature Communications},
+pages = {2831},
+title = {Spatial variation of the rain–snow temperature threshold across the Northern Hemisphere},
+volume = {9},
+year = {2018},
+doi = {10.1038/s41467-018-03629-7},
+}
+
 @article{Jones-Woodward01,
 author = {Jones, J. E. and Woodward, C. S.},
 journal = {Advances in Water Resources},