MODULE soil_therm_inertia !----------------------------------------------------------------------- ! NAME ! soil_therm_inertia ! ! DESCRIPTION ! Compute and update soil thermal properties based on ice content, ! pressure, and cementation state. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, di, k4, minieps ! DECLARATION ! ----------- implicit none ! PARAMETERS ! ---------- real(dp), parameter :: reg_inertie_thresold = 370._dp ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2] real(dp), parameter :: reg_inertie_minvalue = 50._dp ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2] real(dp), parameter :: reg_inertie_maxvalue = 370._dp ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2] real(dp), parameter :: C = 0.0015_dp ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless] real(dp), parameter :: K = 8.1*1e4_dp ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa] real(dp), parameter :: Pa2Torr = 1./133.3_dp ! Conversion from Pa to tor [Pa/Torr] contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE get_ice_TI(ispureice,pore_filling,surf_thermalinertia,ice_thermalinertia) !----------------------------------------------------------------------- ! NAME ! get_ice_TI ! ! DESCRIPTION ! Compute ice thermal properties based on pore filling and purity. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2025 ! ! NOTES ! Uses Siegler et al. (2012) formula for pore-filling ice; ! uses Mellon et al. (2004) value for pure water ice. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use planet, only: porosity ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- logical(k4), intent(in) :: ispureice ! Boolean to check if ice is massive or just pore filling real(dp), intent(in) :: pore_filling ! ice pore filling in each layer (1) real(dp), intent(in) :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2) real(dp), intent(out) :: ice_thermalinertia ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2) ! LOCAL VARIABLES ! --------------- real(dp) :: inertie_purewaterice = 2100._dp ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004) ! CODE ! ---- if (ispureice) then ice_thermalinertia = inertie_purewaterice else ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2) ! Siegler et al., 2012 end if END SUBROUTINE get_ice_TI !======================================================================= !======================================================================= SUBROUTINE update_soil_TI(tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI) !----------------------------------------------------------------------- ! NAME ! update_soil_TI ! ! DESCRIPTION ! Update soil thermal inertia based on ice table, regolith properties, ! and pressure-dependent cementation. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil use soil, only: volcapa, layer, inertiedat, depth_breccia, depth_bedrock, index_breccia, index_bedrock, reg_thprop_dependp use planet, only: TI_breccia, TI_bedrock, TI_regolith_avg, P610 use display, only: print_msg ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: p_avg_new ! Global average surface pressure [Pa] real(dp), dimension(:,:), intent(in) :: tendencies_waterice ! Tendencies on the water ice [kg/m^2/year] real(dp), dimension(:,:), intent(in) :: waterice ! Surface Water ice [kg/m^2] real(dp), dimension(:,:), intent(in) :: icetable_depth ! Ice table depth [m] real(dp), dimension(:,:), intent(in) :: icetable_thickness ! Ice table thickness [m] real(dp), dimension(:,:,:), intent(in) :: ice_porefilling ! Ice porefilling [m^3/m^3] logical(k4), intent(in) :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table real(dp), dimension(:,:,:), intent(inout) :: TI ! Soil Thermal Inertia [J/m^2/K/s^1/2] ! LOCAL VARIABLES ! --------------- integer(di) :: ig, islope, isoil, iref, iend ! Loop variables real(dp), dimension(ngrid,nslope) :: regolith_inertia ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2] real(dp) :: delta ! Difference of depth to compute the thermal inertia in a mixed layer [m] real(dp) :: ice_bottom_depth ! Bottom depth of the subsurface ice [m] real(dp) :: d_part ! Regolith particle size [micrometer] real(dp) :: ice_inertia ! Inertia of water ice [SI] ! CODE ! ---- call print_msg("> Updating soil properties") ! 1. Modification of the regolith thermal inertia. do islope = 1,nslope regolith_inertia(:,islope) = inertiedat(:,1) do ig = 1,ngrid if (tendencies_waterice(ig,islope) < -1.e-5_dp .and. abs(waterice(ig,islope)) < minieps) regolith_inertia(ig,islope) = TI_regolith_avg if (reg_thprop_dependp) then if (TI(ig,1,islope) < reg_inertie_thresold) then d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Torr)**0.6))**(-1./(0.11_dp*log10(P610*Pa2Torr/K))) ! compute particle size, in micrometer regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**0.6*d_part**(-0.11_dp*log10(p_avg_new*Pa2Torr/K))) if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue end if end if end do end do ! 2. Build new Thermal Inertia do islope = 1,nslope do ig = 1,ngrid do isoil = 1,index_breccia TI(ig,isoil,islope) = regolith_inertia(ig,islope) end do if (regolith_inertia(ig,islope) < TI_breccia) then !!! transition delta = depth_breccia TI(ig,index_breccia + 1,islope) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/ & (((delta - layer(index_breccia))/(TI(ig,index_breccia,islope)**2)) + & ((layer(index_breccia + 1) - delta)/(TI_breccia**2)))) do isoil = index_breccia + 2,index_bedrock TI(ig,isoil,islope) = TI_breccia end do else ! we keep the high ti values do isoil = index_breccia + 1,index_bedrock TI(ig,isoil,islope) = TI(ig,index_breccia,islope) end do end if ! TI PEM and breccia comparison !!! transition delta = depth_bedrock TI(ig,index_bedrock + 1,islope) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/ & (((delta - layer(index_bedrock))/(TI(ig,index_bedrock,islope)**2)) + & ((layer(index_bedrock + 1) - delta)/(TI_bedrock**2)))) do isoil = index_bedrock + 2,nsoil TI(ig,isoil,islope) = TI_bedrock end do end do ! ig end do ! islope ! 3. Build new TI in case of ice table do ig = 1,ngrid do islope = 1,nslope if (icetable_depth(ig,islope) > -1.e-6_dp) then ! 3.0 Case where it is perennial ice if (icetable_depth(ig,islope) < minieps) then call get_ice_TI(.true.,1._dp,0._dp,ice_inertia) do isoil = 1,nsoil TI(ig,isoil,islope) = ice_inertia end do else if (icetable_equilibrium) then call get_ice_TI(.false.,1._dp,regolith_inertia(ig,islope),ice_inertia) ! 3.1.1 find the index of the mixed layer iref = 0 ! initialize iref do isoil = 1,nsoil ! loop on layers to find the beginning of the ice table if (icetable_depth(ig,islope) >= layer(isoil)) then iref = isoil ! pure regolith layer up to here else exit ! correct iref was obtained in previous cycle end if end do ! 3.1.2 find the index of the end of the ice table iend = 0 ! initialize iend ice_bottom_depth = icetable_depth(ig,islope) + icetable_thickness(ig,islope) do isoil = 1,nsoil ! loop on layers to find the end of the ice table if (ice_bottom_depth >= layer(isoil)) then iend = isoil ! pure regolith layer up to here else exit ! correct iref was obtained in previous cycle end if end do ! 3.2 Build the new ti if (iref < nsoil) then if (iref == iend) then ! Ice table begins and end in the same mixture with three components if (iref /= 0) then ! mixed layer TI(ig,iref + 1,islope) = sqrt((layer(iref + 1) - layer(iref))/ & (((icetable_depth(ig,islope) - layer(iref))/(TI(ig,iref,islope)**2)) + & ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + & ((layer(iref + 1) - ice_bottom_depth)/(TI(ig,iref + 1,islope)**2)))) else ! first layer is already a mixed layer ! (ie: take layer(iref=0)=0) TI(ig,1,islope) = sqrt((layer(1))/ & (((icetable_depth(ig,islope))/(TI(ig,1,islope)**2)) + & ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + & ((layer(2) - ice_bottom_depth)/(TI(ig,2,islope)**2)))) end if ! of if (iref /= 0) else if (iref /= 0) then ! mixed layer TI(ig,iref + 1,islope) = sqrt((layer(iref + 1) - layer(iref))/ & (((icetable_depth(ig,islope) - layer(iref))/(TI(ig,iref,islope)**2)) + & ((layer(iref + 1) - icetable_depth(ig,islope))/(ice_inertia**2)))) else ! first layer is already a mixed layer ! (ie: take layer(iref=0)=0) TI(ig,1,islope) = sqrt((layer(1))/ & (((icetable_depth(ig,islope))/(TI(ig,1,islope)**2)) + & ((layer(1) - icetable_depth(ig,islope))/(ice_inertia**2)))) end if ! of if (iref /= 0) end if ! iref == iend TI(ig,iref + 2:iend,islope) = ice_inertia if (iend < nsoil) then TI(ig,iend + 1,islope) = sqrt((layer(iend + 1) - layer(iend))/ & (((ice_bottom_depth - layer(iend))/(TI(ig,iend,islope)**2)) + & ((layer(iend + 1) - ice_bottom_depth)/(ice_inertia**2)))) end if end if ! of if (iref < nsoil) else if (icetable_dynamic) then do isoil = 1,nsoil call get_ice_TI(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI(ig,isoil,islope)) end do end if ! of if icetable_equilibrium end if ! permanent glaciers end if ! icetable_depth(ig,islope) > 0. end do !islope end do !ig END SUBROUTINE update_soil_TI !======================================================================= END MODULE soil_therm_inertia