diff --git a/drivers/hrldas/ConfigVarInTransferMod.F90 b/drivers/hrldas/ConfigVarInTransferMod.F90 index f71caa17..7ea82e0f 100644 --- a/drivers/hrldas/ConfigVarInTransferMod.F90 +++ b/drivers/hrldas/ConfigVarInTransferMod.F90 @@ -43,6 +43,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO) noahmp%config%nmlist%OptSurfaceDrag = NoahmpIO%IOPT_SFC noahmp%config%nmlist%OptStomataResistance = NoahmpIO%IOPT_CRS noahmp%config%nmlist%OptSnowAlbedo = NoahmpIO%IOPT_ALB + noahmp%config%nmlist%OptSnowCompact = NoahmpIO%IOPT_COMPACT noahmp%config%nmlist%OptCanopyRadiationTransfer = NoahmpIO%IOPT_RAD noahmp%config%nmlist%OptSnowSoilTempTime = NoahmpIO%IOPT_STC noahmp%config%nmlist%OptSnowThermConduct = NoahmpIO%IOPT_TKSNO @@ -136,7 +137,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO) noahmp%config%domain%VegType = NoahmpIO%NATURAL_TABLE ! set rural vegetation type based on table natural ! urban is handled by explicit urban scheme outside Noah-MP NoahmpIO%GVFMAX(I,J) = 0.96 * 100.0 ! unit: % - endif + endif endif ! treatment for crop point diff --git a/drivers/hrldas/NoahmpIOVarType.F90 b/drivers/hrldas/NoahmpIOVarType.F90 index a936dc04..69719838 100644 --- a/drivers/hrldas/NoahmpIOVarType.F90 +++ b/drivers/hrldas/NoahmpIOVarType.F90 @@ -49,6 +49,7 @@ module NoahmpIOVarType integer :: IOPT_INF ! frozen soil permeability (1-> NY06; 2->Koren99) integer :: IOPT_RAD ! radiation transfer (1->gap=F(3D,cosz); 2->gap=0; 3->gap=1-Fveg) integer :: IOPT_ALB ! snow surface albedo (1->BATS; 2->CLASS) + integer :: IOPT_COMPACT ! snow surface compaction (1->ANDERSON1976; 2->Abolafia-Rosenzweig2024) integer :: IOPT_SNF ! rainfall & snowfall (1-Jordan91; 2->BATS; 3->Noah) integer :: IOPT_TKSNO ! snow thermal conductivity: 1 -> Stieglitz(yen,1965) scheme (default), 2 -> Anderson, 1976 scheme, 3 -> constant, 4 -> Verseghy (1991) scheme, 5 -> Douvill(Yen, 1981) scheme integer :: IOPT_TBOT ! lower boundary of soil temperature (1->zero-flux; 2->Noah) @@ -781,6 +782,12 @@ module NoahmpIOVarType real(kind=kind_noahmp) :: C5_SNOWCOMPACT_TABLE ! snow desctructive metamorphism compaction parameter3 real(kind=kind_noahmp) :: DM_SNOWCOMPACT_TABLE ! upper Limit on destructive metamorphism compaction [kg/m3] real(kind=kind_noahmp) :: ETA0_SNOWCOMPACT_TABLE ! snow viscosity coefficient [kg-s/m2] + real(kind=kind_noahmp) :: SNOWCOMPACTm_TABLE ! snow compaction m parameter for linear sfc temp fitting from AR2024 + real(kind=kind_noahmp) :: SNOWCOMPACTb_TABLE ! snow compaction b parameter for linear sfc temp fitting from AR2024 + real(kind=kind_noahmp) :: SNOWCOMPACT_PSFC1_TABLE ! lower constrain for SnowCompactBurdenFac for high pressure bin from AR2024 + real(kind=kind_noahmp) :: SNOWCOMPACT_PSFC2_TABLE ! lower constrain for SnowCompactBurdenFac for mid pressure bin from AR2024 + real(kind=kind_noahmp) :: SNOWCOMPACT_PSFC3_TABLE ! lower constrain for SnowCompactBurdenFac for low pressure bin from AR2024 + real(kind=kind_noahmp) :: SNOWCOMPACT_Upper_TABLE ! upper constraint on SnowCompactBurdenFac from AR2024 real(kind=kind_noahmp) :: SNLIQMAXFRAC_TABLE ! maximum liquid water fraction in snow real(kind=kind_noahmp) :: SWEMAXGLA_TABLE ! Maximum SWE allowed at glaciers (mm) real(kind=kind_noahmp) :: WSLMAX_TABLE ! maximum lake water storage (mm) diff --git a/drivers/hrldas/NoahmpReadNamelistMod.F90 b/drivers/hrldas/NoahmpReadNamelistMod.F90 index 0b3e1c97..623d253d 100644 --- a/drivers/hrldas/NoahmpReadNamelistMod.F90 +++ b/drivers/hrldas/NoahmpReadNamelistMod.F90 @@ -79,6 +79,7 @@ subroutine NoahmpReadNamelist(NoahmpIO) integer :: frozen_soil_option = 1 integer :: radiative_transfer_option = 3 integer :: snow_albedo_option = 1 + integer :: snow_compaction_option = 2 integer :: snow_thermal_conductivity = 1 integer :: pcp_partition_option = 1 integer :: tbot_option = 2 @@ -123,7 +124,7 @@ subroutine NoahmpReadNamelist(NoahmpIO) forcing_name_LW,forcing_name_SW,forcing_name_PR,forcing_name_SN, & dynamic_veg_option, canopy_stomatal_resistance_option, & btr_option, surface_drag_option, supercooled_water_option, & - frozen_soil_option, radiative_transfer_option, snow_albedo_option, & + frozen_soil_option, radiative_transfer_option, snow_albedo_option, snow_compaction_option, & snow_thermal_conductivity, surface_runoff_option, subsurface_runoff_option, & pcp_partition_option, tbot_option, temp_time_scheme_option, & glacier_option, surface_resistance_option, & @@ -330,6 +331,7 @@ subroutine NoahmpReadNamelist(NoahmpIO) NoahmpIO%IOPT_INF = frozen_soil_option NoahmpIO%IOPT_RAD = radiative_transfer_option NoahmpIO%IOPT_ALB = snow_albedo_option + NoahmpIO%IOPT_COMPACT = snow_compaction_option NoahmpIO%IOPT_SNF = pcp_partition_option NoahmpIO%IOPT_TKSNO = snow_thermal_conductivity NoahmpIO%IOPT_TBOT = tbot_option diff --git a/drivers/hrldas/NoahmpReadTableMod.F90 b/drivers/hrldas/NoahmpReadTableMod.F90 index eb01ceb2..544842a6 100644 --- a/drivers/hrldas/NoahmpReadTableMod.F90 +++ b/drivers/hrldas/NoahmpReadTableMod.F90 @@ -110,14 +110,19 @@ subroutine NoahmpReadTable(NoahmpIO) GRAIN_GROWTH, EXTRA_GROWTH, DIRT_SOOT, BATS_COSZ, BATS_VIS_NEW, & BATS_NIR_NEW, BATS_VIS_AGE, BATS_NIR_AGE, BATS_VIS_DIR, BATS_NIR_DIR, & RSURF_SNOW, RSURF_EXP, C2_SNOWCOMPACT, C3_SNOWCOMPACT, C4_SNOWCOMPACT, & - C5_SNOWCOMPACT, DM_SNOWCOMPACT, ETA0_SNOWCOMPACT, SNLIQMAXFRAC, SWEMAXGLA, & + C5_SNOWCOMPACT, DM_SNOWCOMPACT, ETA0_SNOWCOMPACT, SNOWCOMPACTm,SNOWCOMPACTb, & + SNOWCOMPACT_PSFC1, SNOWCOMPACT_PSFC2, SNOWCOMPACT_PSFC3, SNOWCOMPACT_Upper, & + SNLIQMAXFRAC, SWEMAXGLA, & WSLMAX, ROUS, CMIC, SNOWDEN_MAX, CLASS_ALB_REF, CLASS_SNO_AGE, CLASS_ALB_NEW,& PSIWLT, Z0SOIL, Z0LAKE namelist / noahmp_global_parameters / CO2, O2, TIMEAN, FSATMX, Z0SNO, SSI, SNOW_RET_FAC ,SNOW_EMIS, SWEMX, TAU0, & GRAIN_GROWTH, EXTRA_GROWTH, DIRT_SOOT, BATS_COSZ, BATS_VIS_NEW, & BATS_NIR_NEW, BATS_VIS_AGE, BATS_NIR_AGE, BATS_VIS_DIR, BATS_NIR_DIR, & RSURF_SNOW, RSURF_EXP, C2_SNOWCOMPACT, C3_SNOWCOMPACT, C4_SNOWCOMPACT, & - C5_SNOWCOMPACT, DM_SNOWCOMPACT, ETA0_SNOWCOMPACT, SNLIQMAXFRAC, SWEMAXGLA, & + C5_SNOWCOMPACT, DM_SNOWCOMPACT, ETA0_SNOWCOMPACT, SNOWCOMPACTm,SNOWCOMPACTb, & + SNOWCOMPACT_PSFC1, SNOWCOMPACT_PSFC2, SNOWCOMPACT_PSFC3, SNOWCOMPACT_Upper, & + SNLIQMAXFRAC, SWEMAXGLA, & + SNLIQMAXFRAC, SWEMAXGLA, & WSLMAX, ROUS, CMIC, SNOWDEN_MAX, CLASS_ALB_REF, CLASS_SNO_AGE, CLASS_ALB_NEW,& PSIWLT, Z0SOIL, Z0LAKE @@ -509,6 +514,12 @@ subroutine NoahmpReadTable(NoahmpIO) NoahmpIO%C5_SNOWCOMPACT_TABLE = undefined_real NoahmpIO%DM_SNOWCOMPACT_TABLE = undefined_real NoahmpIO%ETA0_SNOWCOMPACT_TABLE = undefined_real + NoahmpIO%SNOWCOMPACTm_TABLE = undefined_real + NoahmpIO%SNOWCOMPACTb_TABLE = undefined_real + NoahmpIO%SNOWCOMPACT_PSFC1_TABLE = undefined_real + NoahmpIO%SNOWCOMPACT_PSFC2_TABLE = undefined_real + NoahmpIO%SNOWCOMPACT_PSFC3_TABLE = undefined_real + NoahmpIO%SNOWCOMPACT_Upper_TABLE = undefined_real NoahmpIO%SNLIQMAXFRAC_TABLE = undefined_real NoahmpIO%SWEMAXGLA_TABLE = undefined_real NoahmpIO%WSLMAX_TABLE = undefined_real @@ -905,6 +916,12 @@ subroutine NoahmpReadTable(NoahmpIO) NoahmpIO%C5_SNOWCOMPACT_TABLE = C5_SNOWCOMPACT NoahmpIO%DM_SNOWCOMPACT_TABLE = DM_SNOWCOMPACT NoahmpIO%ETA0_SNOWCOMPACT_TABLE = ETA0_SNOWCOMPACT + NoahmpIO%SNOWCOMPACTm_TABLE = SNOWCOMPACTm + NoahmpIO%SNOWCOMPACTb_TABLE = SNOWCOMPACTb + NoahmpIO%SNOWCOMPACT_PSFC1_TABLE = SNOWCOMPACT_PSFC1 + NoahmpIO%SNOWCOMPACT_PSFC2_TABLE = SNOWCOMPACT_PSFC2 + NoahmpIO%SNOWCOMPACT_PSFC3_TABLE = SNOWCOMPACT_PSFC3 + NoahmpIO%SNOWCOMPACT_Upper_TABLE = SNOWCOMPACT_Upper NoahmpIO%SNLIQMAXFRAC_TABLE = SNLIQMAXFRAC NoahmpIO%SWEMAXGLA_TABLE = SWEMAXGLA NoahmpIO%WSLMAX_TABLE = WSLMAX diff --git a/drivers/hrldas/WaterVarInTransferMod.F90 b/drivers/hrldas/WaterVarInTransferMod.F90 index b7208851..0d4f949e 100644 --- a/drivers/hrldas/WaterVarInTransferMod.F90 +++ b/drivers/hrldas/WaterVarInTransferMod.F90 @@ -101,6 +101,12 @@ subroutine WaterVarInTransfer(noahmp, NoahmpIO) noahmp%water%param%SnowCompactAgingFac3 = NoahmpIO%C5_SNOWCOMPACT_TABLE noahmp%water%param%SnowCompactAgingMax = NoahmpIO%DM_SNOWCOMPACT_TABLE noahmp%water%param%SnowViscosityCoeff = NoahmpIO%ETA0_SNOWCOMPACT_TABLE + noahmp%water%param%SnowCompactm = NoahmpIO%SNOWCOMPACTm_TABLE + noahmp%water%param%SnowCompactb = NoahmpIO%SNOWCOMPACTb_TABLE + noahmp%water%param%SnowCompactPSFC1 = NoahmpIO%SNOWCOMPACT_PSFC1_TABLE + noahmp%water%param%SnowCompactPSFC2 = NoahmpIO%SNOWCOMPACT_PSFC2_TABLE + noahmp%water%param%SnowCompactPSFC3 = NoahmpIO%SNOWCOMPACT_PSFC3_TABLE + noahmp%water%param%BurdenFacUpper = NoahmpIO%SNOWCOMPACT_Upper_TABLE noahmp%water%param%SnowLiqFracMax = NoahmpIO%SNLIQMAXFRAC_TABLE noahmp%water%param%SnowLiqHoldCap = NoahmpIO%SSI_TABLE noahmp%water%param%SnowLiqReleaseFac = NoahmpIO%SNOW_RET_FAC_TABLE diff --git a/parameters/NoahmpTable.TBL b/parameters/NoahmpTable.TBL index c9d37c5b..a2246c29 100644 --- a/parameters/NoahmpTable.TBL +++ b/parameters/NoahmpTable.TBL @@ -452,6 +452,12 @@ C5_SNOWCOMPACT = 2.0 ! snow desctructive metamorphism compaction parameter3 DM_SNOWCOMPACT = 100.0 ! upper Limit on destructive metamorphism compaction [kg/m3] ETA0_SNOWCOMPACT = 1.33e+6 ! snow viscosity coefficient [kg-s/m2], Anderson1979: 0.52e6~1.38e6; 1.33e+6 optimized based on SNOTEL obs (He et al. 2021 JGR) + SNOWCOMPACTm = -0.00069503 ! snow compaction m parameter for linear sfc temp fitting from AR2024 + SNOWCOMPACTb = 0.20606699 ! snow compaction b parameter for linear sfc temp fitting from AR2024 + SNOWCOMPACT_PSFC1 = 0.017 ! lower constrain for SnowCompactBurdenFac for high pressure bin from AR2024 + SNOWCOMPACT_PSFC2 = 0.018 ! lower constrain for SnowCompactBurdenFac for mid pressure bin from AR2024 + SNOWCOMPACT_PSFC3 = 0.019 ! lower constrain for SnowCompactBurdenFac for low pressure bin from AR2024 + SNOWCOMPACT_Upper = 0.0315 ! upper constraint on SnowCompactBurdenFac from AR2024 SNLIQMAXFRAC = 0.4 ! maximum liquid water fraction in snow SWEMAXGLA = 5000.0 ! Maximum SWE allowed at glaciers (mm) SNOWDEN_MAX = 120.0 ! maximum fresh snowfall density (kg/m3) @@ -460,6 +466,7 @@ CLASS_ALB_NEW = 0.84 ! fresh snow albedo in CLASS scheme RSURF_SNOW = 50.0 ! surface resistence for snow [s/m] Z0SNO = 0.002 ! snow surface roughness length (m) + ! other soil and hydrological parameters RSURF_EXP = 5.0 ! exponent in the shape parameter for soil resistance option 1 WSLMAX = 5000.0 ! maximum lake water storage (mm) diff --git a/src/ConfigVarInitMod.F90 b/src/ConfigVarInitMod.F90 index 5c8af537..68d7f927 100644 --- a/src/ConfigVarInitMod.F90 +++ b/src/ConfigVarInitMod.F90 @@ -30,6 +30,7 @@ subroutine ConfigVarInitDefault(noahmp) noahmp%config%nmlist%OptSurfaceDrag = undefined_int noahmp%config%nmlist%OptStomataResistance = undefined_int noahmp%config%nmlist%OptSnowAlbedo = undefined_int + noahmp%config%nmlist%OptSnowCompact = undefined_int noahmp%config%nmlist%OptCanopyRadiationTransfer = undefined_int noahmp%config%nmlist%OptSnowSoilTempTime = undefined_int noahmp%config%nmlist%OptSnowThermConduct = undefined_int diff --git a/src/ConfigVarType.F90 b/src/ConfigVarType.F90 index dc7979f3..0b59d887 100644 --- a/src/ConfigVarType.F90 +++ b/src/ConfigVarType.F90 @@ -121,6 +121,10 @@ module ConfigVarType ! 1 -> include phase change of ice (default) ! 2 -> ice treatment more like original Noah + integer :: OptSnowCompact ! options for ground snow compaction + ! 1 -> original scheme from Anderson (1976) + ! 2 -> enhanced scheme from Abolafia-Rosenzweig et al. (2024) + end type namelist_type diff --git a/src/Makefile b/src/Makefile index 706c0219..eca3f701 100644 --- a/src/Makefile +++ b/src/Makefile @@ -42,6 +42,7 @@ OBJS = ConstantDefineMod.o \ SnowLayerDivideMod.o \ SnowLayerWaterComboMod.o \ SnowpackCompactionMod.o \ + SnowpackCompactionARMod.o \ SnowpackHydrologyMod.o \ SnowWaterMainMod.o \ SoilHydraulicPropertyMod.o \ @@ -205,10 +206,11 @@ SnowLayerDivideMod.o: ../utility/Machine.o NoahmpVarType.o Const SnowLayerWaterComboMod.o SnowLayerWaterComboMod.o: ../utility/Machine.o ConstantDefineMod.o SnowpackCompactionMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o +SnowpackCompactionARMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SnowpackHydrologyMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o \ SnowLayerCombineMod.o SnowWaterMainMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o \ - SnowfallBelowCanopyMod.o SnowpackCompactionMod.o SnowLayerDivideMod.o \ + SnowfallBelowCanopyMod.o SnowpackCompactionMod.o SnowpackCompactionARMod.o SnowLayerDivideMod.o \ SnowLayerCombineMod.o SnowpackHydrologyMod.o SoilHydraulicPropertyMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SoilMoistureSolverMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o \ @@ -339,7 +341,7 @@ PsychrometricVariableGlacierMod.o: ../utility/Machine.o NoahmpVarType.o Co ResistanceGroundEvaporationGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SnowCoverGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SnowWaterMainGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SnowfallBelowCanopyMod.o \ - SnowpackCompactionMod.o SnowLayerCombineMod.o SnowLayerDivideMod.o \ + SnowpackCompactionMod.o SnowpackCompactionARMod.o SnowLayerCombineMod.o SnowLayerDivideMod.o \ SnowpackHydrologyGlacierMod.o SnowpackHydrologyGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SnowLayerCombineMod.o SurfaceAlbedoGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SnowAgingBatsMod.o \ diff --git a/src/SnowWaterMainMod.F90 b/src/SnowWaterMainMod.F90 index 2e3e7f00..24d89be2 100644 --- a/src/SnowWaterMainMod.F90 +++ b/src/SnowWaterMainMod.F90 @@ -8,6 +8,7 @@ module SnowWaterMainMod use ConstantDefineMod use SnowfallBelowCanopyMod, only : SnowfallAfterCanopyIntercept use SnowpackCompactionMod, only : SnowpackCompaction + use SnowpackCompactionARMod, only : SnowpackCompactionAR use SnowLayerCombineMod, only : SnowLayerCombine use SnowLayerDivideMod, only : SnowLayerDivide use SnowpackHydrologyMod, only : SnowpackHydrology @@ -39,6 +40,7 @@ subroutine SnowWaterMain(noahmp) MainTimeStep => noahmp%config%domain%MainTimeStep ,& ! in, noahmp main time step [s] DepthSoilLayer => noahmp%config%domain%DepthSoilLayer ,& ! in, depth [m] of layer-bottom from soil surface SnoWatEqvMaxGlacier => noahmp%water%param%SnoWatEqvMaxGlacier ,& ! in, Maximum SWE allowed at glaciers [mm] + OptSnowCompact => noahmp%config%nmlist%OptSnowCompact ,& ! in, options for ground snow compaction ThicknessSnowSoilLayer => noahmp%config%domain%ThicknessSnowSoilLayer ,& ! inout, thickness of snow/soil layers [m] DepthSnowSoilLayer => noahmp%config%domain%DepthSnowSoilLayer ,& ! inout, depth of snow/soil layer-bottom [m] NumSnowLayerNeg => noahmp%config%domain%NumSnowLayerNeg ,& ! inout, actual number of snow layers (negative) @@ -64,7 +66,8 @@ subroutine SnowWaterMain(noahmp) ! do following snow layer compaction, combination, and division only for multi-layer snowpack ! snowpack compaction - if ( NumSnowLayerNeg < 0 ) call SnowpackCompaction(noahmp) + if ( NumSnowLayerNeg < 0 .and. OptSnowCompact == 1) call SnowpackCompaction(noahmp) + if ( NumSnowLayerNeg < 0 .and. OptSnowCompact == 2) call SnowpackCompactionAR(noahmp) ! snow layer combination if ( NumSnowLayerNeg < 0 ) call SnowLayerCombine(noahmp) diff --git a/src/SnowpackCompactionARMod.F90 b/src/SnowpackCompactionARMod.F90 new file mode 100644 index 00000000..0e104790 --- /dev/null +++ b/src/SnowpackCompactionARMod.F90 @@ -0,0 +1,150 @@ +module SnowpackCompactionARMod + +!!! Snowpack compaction process +!!! Update snow depth via compaction due to destructive metamorphism, overburden, & melt + + use Machine + use NoahmpVarType + use ConstantDefineMod + + implicit none + +contains + + subroutine SnowpackCompactionAR(noahmp) + +! ------------------------ Code history ----------------------------------- +! Original Noah-MP subroutine: COMPACT +! Original code: Guo-Yue Niu and Noah-MP team (Niu et al. 2011) +! Refactered code: C. He, P. Valayamkunnath, & refactor team (He et al. 2023) +! Enhanced snow compaction scheme: R. Abolafia-Rosenzweig (Abolafia-Rosenzweig et al. 2024) +! ------------------------------------------------------------------------- + + implicit none + + type(noahmp_type), intent(inout) :: noahmp + +! local variable + integer :: LoopInd ! snow layer loop index + real(kind=kind_noahmp) :: SnowBurden ! pressure of overlying snow [kg/m2] + real(kind=kind_noahmp) :: SnowCompactAgeExpFac ! EXPF=exp(-c4*(273.15-TemperatureSoilSnow)) + real(kind=kind_noahmp) :: TempDiff ! ConstFreezePoint - TemperatureSoilSnow[K] + real(kind=kind_noahmp) :: SnowVoid ! void (1 - SnowIce - SnowLiqWater) + real(kind=kind_noahmp) :: SnowWatTotTmp ! water mass (ice + liquid) [kg/m2] + real(kind=kind_noahmp) :: SnowIceDens ! partial density of ice [kg/m3] + +! -------------------------------------------------------------------- + associate( & + MainTimeStep => noahmp%config%domain%MainTimeStep ,& ! in, noahmp main time step [s] + TemperatureSoilSnow => noahmp%energy%state%TemperatureSoilSnow ,& ! in, snow and soil layer temperature [K] + SnowIce => noahmp%water%state%SnowIce ,& ! in, snow layer ice [mm] + SnowLiqWater => noahmp%water%state%SnowLiqWater ,& ! in, snow layer liquid water [mm] + IndexPhaseChange => noahmp%water%state%IndexPhaseChange ,& ! in, phase change index [0-none;1-melt;2-refreeze] + SnowIceFracPrev => noahmp%water%state%SnowIceFracPrev ,& ! in, ice fraction in snow layers at previous timestep + SnowCompactBurdenFac => noahmp%water%param%SnowCompactBurdenFac ,& ! in, snow overburden compaction parameter [m3/kg] + SnowCompactAgingFac1 => noahmp%water%param%SnowCompactAgingFac1 ,& ! in, snow desctructive metamorphism compaction factor1 [1/s] + SnowCompactAgingFac2 => noahmp%water%param%SnowCompactAgingFac2 ,& ! in, snow desctructive metamorphism compaction factor2 [1/k] + SnowCompactAgingFac3 => noahmp%water%param%SnowCompactAgingFac3 ,& ! in, snow desctructive metamorphism compaction factor3 + SnowCompactm => noahmp%water%param%SnowCompactm ,& ! snow compaction m parameter for linear sfc temp fitting from AR2024 + SnowCompactb => noahmp%water%param%SnowCompactb ,& ! snow compaction b parameter for linear sfc temp fitting from AR2024 + SnowCompactPSFC1 => noahmp%water%param%SnowCompactPSFC1 ,& ! lower constrain for SnowCompactBurdenFac for high pressure bin from AR2024 + SnowCompactPSFC2 => noahmp%water%param%SnowCompactPSFC2 ,& ! lower constrain for SnowCompactBurdenFac for mid pressure bin from AR2024 + SnowCompactPSFC3 => noahmp%water%param%SnowCompactPSFC3 ,& ! lower constrain for SnowCompactBurdenFac for low pressure bin from AR2024 + BurdenFacUpper => noahmp%water%param%BurdenFacUpper ,& ! upper constraint on SnowCompactBurdenFac from AR2024 + SnowCompactAgingMax => noahmp%water%param%SnowCompactAgingMax ,& ! in, maximum destructive metamorphism compaction [kg/m3] + SnowViscosityCoeff => noahmp%water%param%SnowViscosityCoeff ,& ! in, snow viscosity coeff [kg s/m2],Anderson1979:0.52e6~1.38e6 + NumSnowLayerNeg => noahmp%config%domain%NumSnowLayerNeg ,& ! inout, actual number of snow layers (negative) + ThicknessSnowSoilLayer => noahmp%config%domain%ThicknessSnowSoilLayer ,& ! inout, thickness of snow/soil layers [m] + CompactionSnowAging => noahmp%water%flux%CompactionSnowAging ,& ! out, rate of compaction due to destructive metamorphism [1/s] + CompactionSnowBurden => noahmp%water%flux%CompactionSnowBurden ,& ! out, rate of compaction of snowpack due to overburden [1/s] + CompactionSnowMelt => noahmp%water%flux%CompactionSnowMelt ,& ! out, rate of compaction of snowpack due to melt [1/s] + CompactionSnowTot => noahmp%water%flux%CompactionSnowTot ,& ! out, change in fractional-thickness due to compaction [1/s] + SnowIceFrac => noahmp%water%state%SnowIceFrac ,& ! out, fraction of ice in snow layers at current time step + TemperatureAirRefHeight => noahmp%forcing%TemperatureAirRefHeight ,& ! in, air temperature [K] at reference height + PressureAirRefHeight => noahmp%forcing%PressureAirRefHeight & ! in, air pressure [Pa] at reference height + ) +! ---------------------------------------------------------------------- + +! initialization for out-only variables + CompactionSnowAging(:) = 0.0 + CompactionSnowBurden(:) = 0.0 + CompactionSnowMelt(:) = 0.0 + CompactionSnowTot(:) = 0.0 + SnowIceFrac(:) = 0.0 + + ! start snow compaction + SnowBurden = 0.0 + !SnowCompactBurdenFac updated from Abolafia-Rosenzweig et al., 2024 + SnowCompactBurdenFac = SnowCompactm * TemperatureAirRefHeight + SnowCompactb + !pressure-based lower constraints: + if (PressureAirRefHeight>=85000) then !high pressure bin + SnowCompactBurdenFac = max(SnowCompactBurdenFac,SnowCompactPSFC1) + endif + if (PressureAirRefHeight>=80000 .and. PressureAirRefHeight<85000) then !mid-pressure bin + SnowCompactBurdenFac = max(SnowCompactBurdenFac,SnowCompactPSFC2) + endif + if (PressureAirRefHeight<80000) then !low pressure bin + SnowCompactBurdenFac = max(SnowCompactBurdenFac,SnowCompactPSFC3) + endif + !upper constraint on SnowCompactBurdenFac + SnowCompactBurdenFac = min(SnowCompactBurdenFac,BurdenFacUpper) + + do LoopInd = NumSnowLayerNeg+1, 0 + + SnowWatTotTmp = SnowIce(LoopInd) + SnowLiqWater(LoopInd) + SnowIceFrac(LoopInd) = SnowIce(LoopInd) / SnowWatTotTmp + SnowVoid = 1.0 - (SnowIce(LoopInd)/ConstDensityIce + SnowLiqWater(LoopInd)/ConstDensityWater) / & + ThicknessSnowSoilLayer(LoopInd) + + ! Allow compaction only for non-saturated node and higher ice lens node. + if ( (SnowVoid > 0.001) .and. (SnowIce(LoopInd) > 0.1) ) then + SnowIceDens = SnowIce(LoopInd) / ThicknessSnowSoilLayer(LoopInd) + TempDiff = max(0.0, ConstFreezePoint-TemperatureSoilSnow(LoopInd)) + + ! Settling/compaction as a result of destructive metamorphism + SnowCompactAgeExpFac = exp(-SnowCompactAgingFac2 * TempDiff) + CompactionSnowAging(LoopInd) = -SnowCompactAgingFac1 * SnowCompactAgeExpFac + if ( SnowIceDens > SnowCompactAgingMax ) & + CompactionSnowAging(LoopInd) = CompactionSnowAging(LoopInd) * exp(-46.0e-3*(SnowIceDens-SnowCompactAgingMax)) + if ( SnowLiqWater(LoopInd) > (0.01*ThicknessSnowSoilLayer(LoopInd)) ) & + CompactionSnowAging(LoopInd) = CompactionSnowAging(LoopInd) * SnowCompactAgingFac3 ! Liquid water term + + ! Compaction due to overburden + CompactionSnowBurden(LoopInd) = -(SnowBurden + 0.5*SnowWatTotTmp) * & + exp(-0.08*TempDiff-SnowCompactBurdenFac*SnowIceDens) / SnowViscosityCoeff ! 0.5*SnowWatTotTmp -> self-burden + + ! Compaction occurring during melt + if ( IndexPhaseChange(LoopInd) == 1 ) then + CompactionSnowMelt(LoopInd) = max(0.0, (SnowIceFracPrev(LoopInd)-SnowIceFrac(LoopInd)) / & + max(1.0e-6, SnowIceFracPrev(LoopInd))) + CompactionSnowMelt(LoopInd) = -CompactionSnowMelt(LoopInd) / MainTimeStep ! sometimes too large + else + CompactionSnowMelt(LoopInd) = 0.0 + endif + + ! Time rate of fractional change in snow thickness (units of s-1) + CompactionSnowTot(LoopInd) = (CompactionSnowAging(LoopInd) + CompactionSnowBurden(LoopInd) + & + CompactionSnowMelt(LoopInd) ) * MainTimeStep + CompactionSnowTot(LoopInd) = max(-0.5, CompactionSnowTot(LoopInd)) + + ! The change in DZ due to compaction + ThicknessSnowSoilLayer(LoopInd) = ThicknessSnowSoilLayer(LoopInd) * (1.0 + CompactionSnowTot(LoopInd)) + ThicknessSnowSoilLayer(LoopInd) = max(ThicknessSnowSoilLayer(LoopInd), & + SnowIce(LoopInd)/ConstDensityIce + SnowLiqWater(LoopInd)/ConstDensityWater) + + ! Constrain snow density to a reasonable range (50~500 kg/m3) + ThicknessSnowSoilLayer(LoopInd) = min( max( ThicknessSnowSoilLayer(LoopInd),& + (SnowIce(LoopInd)+SnowLiqWater(LoopInd))/500.0 ), & + (SnowIce(LoopInd)+SnowLiqWater(LoopInd))/50.0 ) + endif + + ! Pressure of overlying snow + SnowBurden = SnowBurden + SnowWatTotTmp + + enddo + + end associate + + end subroutine SnowpackCompactionAR + +end module SnowpackCompactionARMod diff --git a/src/WaterVarInitMod.F90 b/src/WaterVarInitMod.F90 index a03d8b4f..22de1b11 100644 --- a/src/WaterVarInitMod.F90 +++ b/src/WaterVarInitMod.F90 @@ -230,6 +230,12 @@ subroutine WaterVarInitDefault(noahmp) noahmp%water%param%SnowCompactAgingFac3 = undefined_real noahmp%water%param%SnowCompactAgingMax = undefined_real noahmp%water%param%SnowViscosityCoeff = undefined_real + noahmp%water%param%SnowCompactm = undefined_real + noahmp%water%param%SnowCompactb = undefined_real + noahmp%water%param%SnowCompactPSFC1 = undefined_real + noahmp%water%param%SnowCompactPSFC2 = undefined_real + noahmp%water%param%SnowCompactPSFC3 = undefined_real + noahmp%water%param%BurdenFacUpper = undefined_real noahmp%water%param%SnowLiqFracMax = undefined_real noahmp%water%param%SnowLiqHoldCap = undefined_real noahmp%water%param%SnowLiqReleaseFac = undefined_real diff --git a/src/WaterVarType.F90 b/src/WaterVarType.F90 index 2d2f9132..39443521 100644 --- a/src/WaterVarType.F90 +++ b/src/WaterVarType.F90 @@ -173,6 +173,12 @@ module WaterVarType real(kind=kind_noahmp) :: SnowCompactAgingFac3 ! snow desctructive metamorphism compaction parameter3 real(kind=kind_noahmp) :: SnowCompactAgingMax ! upper Limit on destructive metamorphism compaction [kg/m3] real(kind=kind_noahmp) :: SnowViscosityCoeff ! snow viscosity coefficient [kg-s/m2], Anderson1979: 0.52e6~1.38e6 + real(kind=kind_noahmp) :: SnowCompactm ! snow compaction m parameter for linear sfc temp fitting from AR2024 + real(kind=kind_noahmp) :: SnowCompactb ! snow compaction b parameter for linear sfc temp fitting from AR2024 + real(kind=kind_noahmp) :: SnowCompactPSFC1 ! lower constrain for SnowCompactBurdenFac for high pressure bin from AR2024 + real(kind=kind_noahmp) :: SnowCompactPSFC2 ! lower constrain for SnowCompactBurdenFac for mid pressure bin from AR2024 + real(kind=kind_noahmp) :: SnowCompactPSFC3 ! lower constrain for SnowCompactBurdenFac for low pressure bin from AR2024 + real(kind=kind_noahmp) :: BurdenFacUpper ! upper constraint on SnowCompactBurdenFac from AR2024 real(kind=kind_noahmp) :: SnowLiqFracMax ! maximum liquid water fraction in snow real(kind=kind_noahmp) :: SnowLiqHoldCap ! liquid water holding capacity for snowpack [m3/m3] real(kind=kind_noahmp) :: SnowLiqReleaseFac ! snowpack water release timescale factor [1/s]