MODULE sorption !----------------------------------------------------------------------- ! NAME ! sorption ! ! DESCRIPTION ! CO2 and H2O adsorption/desorption in regolith following Zent & Quinn ! (1995) and Jackosky et al. (1997). ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 02/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, qp, di, k4, minieps ! DECLARATION ! ----------- implicit none ! PARAMETERS ! ---------- logical(k4), protected :: do_sorption ! Flag to compute adsorption/desorption contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE set_sorption_config(do_sorption_in) !----------------------------------------------------------------------- ! NAME ! set_sorption_config ! ! DESCRIPTION ! Setter for 'sorption' configuration parameters. ! ! AUTHORS & DATE ! JB Clement, 02/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use utility, only: bool2str use display, only: print_msg ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- logical(k4), intent(in) :: do_sorption_in ! CODE ! ---- do_sorption = do_sorption_in call print_msg('do_sorption = '//bool2str(do_sorption)) END SUBROUTINE set_sorption_config !======================================================================= !======================================================================= SUBROUTINE evolve_regolith_adsorption(d_h2oice,d_co2ice,h2oice,co2ice,tsoil,TI,ps,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads) !----------------------------------------------------------------------- ! NAME ! evolve_regolith_adsorption ! ! DESCRIPTION ! Compute both H2O and CO2 adsorption in regolith. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 12/2025 ! C. Metz, 02/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(in) :: d_h2oice, d_co2ice ! tendencies on the glaciers [1] real(dp), dimension(:,:), intent(in) :: h2oice ! H2O ice at the surface [kg/m^2] real(dp), dimension(:,:), intent(in) :: co2ice ! CO2 ice at the surface [kg/m^2] real(dp), dimension(:,:,:), intent(in) :: TI ! Soil Thermal inertia [J/K/^2/s^1/2] real(dp), dimension(:,:,:), intent(in) :: tsoil ! Soil temperature [K] real(dp), dimension(:,:), intent(in) :: ps ! Average surface pressure [Pa] real(dp), dimension(:,:), intent(in) :: q_co2_ts ! Mass mixing ratio of CO2 in the first layer [kg/kg] real(dp), dimension(:,:), intent(in) :: q_h2o_ts ! Mass mixing ratio of H2O in the first layer [kg/kg] real(dp), dimension(:), intent(out) :: delta_h2o_ads ! Difference density of H2O adsorbed [kg/m^3] real(dp), dimension(:), intent(out) :: delta_co2_ads ! Difference density of CO2 adsorbed [kg/m^3] real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg, co2_ads_reg ! LOCAL VARIABLES ! --------------- real(dp), dimension(ngrid,nsoil,nslope) :: theta_h2o_ads ! Fraction of the pores occupied by H2O molecules ! CODE ! ---- call evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,theta_h2o_ads,h2o_ads_reg,delta_h2o_ads) call evolve_co2_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads) END SUBROUTINE evolve_regolith_adsorption !======================================================================= !======================================================================= SUBROUTINE evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,theta_h2o_ads,h2o_ads_reg,delta_h2o_ads) !----------------------------------------------------------------------- ! NAME ! evolve_h2o_ads ! ! DESCRIPTION ! Compute H2O adsorption in regolith following Jackosky et al. (1997). ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil, nday use soil, only: layer, index_breccia use slopes, only: subslope_dist, def_slope_mean use atmosphere, only: ap, bp use planet, only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2,m_noco2, rho_regolith use maths, only: pi ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(in) :: d_h2oice, d_co2ice ! Tendencies on the glaciers () real(dp), dimension(:,:), intent(in) :: h2oice ! H2O ice at the surface [kg/m^2] real(dp), dimension(:,:), intent(in) :: co2ice ! CO2 ice at the surface [kg/m^2] real(dp), dimension(:,:), intent(in) :: ps ! Surface pressure [Pa] real(dp), dimension(:,:), intent(in) :: q_co2_ts ! Mass mixing ratio of CO2 in the first layer [kg/kg] real(dp), dimension(:,:), intent(in) :: q_h2o_ts ! Mass mixing ratio of H2O in the first layer [kg/kg] real(dp), dimension(:,:,:), intent(in) :: TI ! Soil Thermal inertia [J/K/^2/s^1/2] real(dp), dimension(:,:,:), intent(in) :: tsoil ! Soil temperature [K] real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg ! Density of H2O adsorbed [kg/m^3] real(dp), dimension(:,:,:), intent(out) :: theta_h2o_ads ! Fraction of the pores occupied by H2O molecules real(dp), dimension(:), intent(out) :: delta_h2o_ads ! Difference density of H2O adsorbed [kg/m^3] ! LOCAL VARIABLES ! --------------- ! Constants: real(dp) :: Ko = 1.57e-8_dp ! Jackosky et al. 1997 real(dp) :: e = 2573.9_dp ! Jackosky et al. 1997 real(dp) :: mu = 0.48_dp ! Jackosky et al. 1997 real(dp) :: m_theta = 2.84e-7_dp ! Mass of H2O per m^2 absorbed Jackosky et al. 1997 ! real(dp) :: as = 18.9e3_dp ! Specific area, Buhler & Piqueux 2021 real(dp) :: as = 9.48e4_dp ! Specific area, Zent real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation real(dp), dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete ! Difference in the mass per slope and soil layer [kg/m^3] real(dp) :: K ! Used to compute theta integer(di) :: ig, iloop, islope, it ! For loops logical(k4), dimension(ngrid,nslope) :: ispermanent_co2glaciers ! Check if the CO2 glacier is permanent logical(k4), dimension(ngrid,nslope) :: ispermanent_h2oglaciers ! Check if the H2O glacier is permanent real(dp), dimension(ngrid,nslope) :: deltam_reg_slope ! Difference density of H2O adsorbed per slope [kg/m^3] real(dp), dimension(ngrid,nsoil,nslope) :: dm_h2o_regolith_slope ! Elementary H2O mass adsorded per mesh per slope real(dp) :: A, B ! Used to compute the mean mass above the surface real(dp) :: p_sat ! Saturated vapor pressure of ice real(dp), dimension(:,:), allocatable :: mass_mean ! Mean mass above the surface real(dp), dimension(:,:), allocatable :: zplev_mean ! Pressure above the surface real(dp), dimension(:,:), allocatable :: pvapor ! Partial pressure above the surface real(dp), dimension(:) , allocatable :: pvapor_avg ! Yearly averaged ! CODE ! ---- ! 0. Some initializations allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pvapor(ngrid,nday),pvapor_avg(ngrid)) A = 1._dp/m_co2 - 1._dp/m_noco2 B = 1._dp/m_noco2 theta_h2o_ads = 0._dp dm_h2o_regolith_slope = 0._dp ispermanent_h2oglaciers = .false. ispermanent_co2glaciers = .false. ! 0.1 Look at perennial ice where (abs(d_h2oice(:,:)) > 1.e-5 .and. h2oice(:,:) > 0._dp) ispermanent_h2oglaciers(:,:) = .true. where (abs(d_co2ice(:,:)) > 1.e-5 .and. co2ice(:,:) > 0._dp) ispermanent_co2glaciers(:,:) = .true. ! 0.2 Compute the partial pressure of vapor ! a. the molecular mass into the column do ig = 1,ngrid mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B) end do ! b. pressure level do it = 1,nday do ig = 1,ngrid zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it) end do end do ! c. Vapor pressure pvapor = mass_mean/m_h2o*q_h2o_ts*zplev_mean pvapor_avg = sum(pvapor,2)/nday deallocate(pvapor,zplev_mean,mass_mean) ! 1. we compute the mass of H2O adsorded in each layer of the meshes do ig = 1,ngrid do islope = 1,nslope do iloop = 1,index_breccia K = Ko*exp(e/tsoil(ig,iloop,islope)) if (TI(ig,iloop,islope) < inertie_thresold) then theta_h2o_ads(ig,iloop,islope) = (K*pvapor_avg(ig)/(1._dp + K*pvapor_avg(ig)))**mu else p_sat = exp(beta_clap_h2o/tsoil(ig,iloop,islope) + alpha_clap_h2o) ! we assume fixed temperature in the ice ... not really good but ... theta_h2o_ads(ig,iloop,islope) = (K*p_sat/(1._dp + K*p_sat))**mu end if dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_ads(ig,iloop,islope)*m_theta*rho_regolith end do end do end do ! 2. Check the exchange between the atmosphere and the regolith do ig = 1,ngrid delta_h2o_ads(ig) = 0._dp do islope = 1,nslope deltam_reg_slope(ig,islope) = 0._dp do iloop = 1,index_breccia if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then if (iloop == 1) then deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop)) else deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1)) end if else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC! deltam_reg_complete(ig,iloop,islope) = 0._dp end if deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope) end do delta_h2o_ads(ig) = delta_h2o_ads(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180._dp) end do end do h2o_ads_reg = dm_h2o_regolith_slope END SUBROUTINE evolve_h2o_ads !======================================================================= !======================================================================= SUBROUTINE evolve_co2_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads) !----------------------------------------------------------------------- ! NAME ! evolve_co2_ads ! ! DESCRIPTION ! Compute CO2 adsorption in regolith following Zent & Quinn (1995). ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil, nday use soil, only: layer, index_breccia use slopes, only: subslope_dist, def_slope_mean use atmosphere, only: ap, bp use planet, only: m_co2, m_noco2, rho_regolith use maths, only: pi ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(in) :: ps ! Average surface pressure [Pa] real(dp), dimension(:,:,:), intent(in) :: tsoil ! Mean Soil Temperature [K] real(dp), dimension(:,:,:), intent(in) :: TI ! Mean Thermal Inertia [USI] real(dp), dimension(:,:), intent(in) :: d_h2oice, d_co2ice ! Tendencies on the glaciers () real(dp), dimension(:,:), intent(in) :: q_co2_ts, q_h2o_ts ! Mass mixing ratio of CO2 and H2O in the first layer [kg/kg] real(dp), dimension(:,:), intent(in) :: h2oice ! Water ice at the surface [kg/m^2] real(dp), dimension(:,:), intent(in) :: co2ice ! CO2 ice at the surface [kg/m^2] real(dp), dimension(:,:,:), intent(inout) :: co2_ads_reg ! Density of CO2 adsorbed [kg/m^3] real(dp), dimension(:), intent(out) :: delta_co2_ads ! Difference density of CO2 adsorbed [kg/m^3] ! LOCAL VARIABLES ! --------------- ! Constants: real(dp) :: alpha = 7.512e-6_dp ! Zent & Quinn 1995 real(dp) :: beta = -1541.5_dp ! Zent & Quinn 1995 real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation real(dp) :: m_theta = 4.27e-7_dp ! Mass of CO2 per m^2 absorbed ! real(dp) :: as = 18.9e3_dp ! Specific area, Buhler & Piqueux 2021 real(dp) :: as = 9.48e4_dp ! Same as previous but from zent real(dp) :: A, B ! Used to compute the mean mass above the surface integer(di) :: ig, islope, iloop, it ! For loops real(dp), dimension(ngrid,nsoil,nslope) :: dm_co2_regolith_slope ! Elementary mass adsorded per mesh per slope logical(k4), dimension(ngrid,nslope) :: ispermanent_co2glaciers ! Check if the CO2 glacier is permanent logical(k4), dimension(ngrid,nslope) :: ispermanent_h2oglaciers ! Check if the H2O glacier is permanent real(dp), dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete ! Difference in the mass per slope and soil layer [kg/m^3] real(dp), dimension(ngrid,nslope) :: deltam_reg_slope ! Difference in the mass per slope [kg/m^3] real(dp), dimension(ngrid,nsoil,nslope) :: m_h2o_adsorbed ! Density of CO2 adsorbed [kg/m^3] real(dp), dimension(ngrid,nsoil,nslope) :: theta_h2o_ads ! Fraction of the pores occupied by H2O molecules real(dp), dimension(ngrid) :: delta_mh2o ! Difference density of H2O adsorbed [kg/m^3] ! nday array are allocated because heavy... real(dp), dimension(:,:), allocatable :: mass_mean ! Mean mass above the surface real(dp), dimension(:,:), allocatable :: zplev_mean ! Pressure above the surface real(dp), dimension(:,:), allocatable :: pco2 ! Partial pressure above the surface real(dp), dimension(:), allocatable :: pco2_avg ! Yearly averaged ! CODE ! ---- ! 0. Some initializations allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pco2(ngrid,nday),pco2_avg(ngrid)) m_h2o_adsorbed = 0._dp A = 1._dp/m_co2 - 1._dp/m_noco2 B = 1._dp/m_noco2 dm_co2_regolith_slope = 0._dp delta_co2_ads = 0._dp ispermanent_h2oglaciers = .false. ispermanent_co2glaciers = .false. ! 0.1 Look at perennial ice where (abs(d_h2oice(:,:)) > 1.e-5 .and. h2oice(:,:) > 0._dp) ispermanent_h2oglaciers(:,:) = .true. where (abs(d_co2ice(:,:)) > 1.e-5 .and. co2ice(:,:) > 0._dp) ispermanent_co2glaciers(:,:) = .true. ! 0.2 Compute the partial pressure of CO2 ! a. the molecular mass into the column do ig = 1,ngrid mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B) end do ! b. pressure level do it = 1,nday do ig = 1,ngrid zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it) end do end do ! c. Vapor pressure pco2 = mass_mean/m_co2*q_co2_ts*zplev_mean pco2_avg(:) = sum(pco2(:,:),2)/nday deallocate(zplev_mean,mass_mean,pco2) ! 1. Compute the fraction of the pores occupied by H2O call evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_co2_ts,q_h2o_ts,tsoil,TI,theta_h2o_ads,m_h2o_adsorbed,delta_mh2o) ! 2. we compute the mass of CO2 adsorded in each layer of the meshes do ig = 1,ngrid do islope = 1,nslope do iloop = 1,index_breccia if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig)/ & (alpha*pco2_avg(ig) + sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope))) else if (abs(co2_ads_reg(ig,iloop,islope)) < minieps) then !!! we are at first call dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig) & /(alpha*pco2_avg(ig)+sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope))) else ! no change: permanent ice stick the atoms of CO2 dm_co2_regolith_slope(ig,iloop,islope) = co2_ads_reg(ig,iloop,islope) end if end if end do end do end do ! 3. Check the exchange between the atmosphere and the regolith do ig = 1,ngrid delta_co2_ads(ig) = 0._dp do islope = 1,nslope deltam_reg_slope(ig,islope) = 0._dp do iloop = 1,index_breccia if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then if (iloop == 1) then deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop)) else deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1)) end if else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC! deltam_reg_complete(ig,iloop,islope) = 0._dp end if deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope) end do delta_co2_ads(ig) = delta_co2_ads(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180._dp) end do end do co2_ads_reg = dm_co2_regolith_slope END SUBROUTINE evolve_co2_ads !======================================================================= !======================================================================= SUBROUTINE compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2,totmass_adsh2o) !----------------------------------------------------------------------- ! NAME ! compute_totmass_adsorbed ! ! DESCRIPTION ! Compute the total mass of CO2 and H2O adsorbed in the regolith. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 12/2025 ! C. Metz, 02/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil, cell_area use slopes, only: subslope_dist, def_slope_mean use soil, only: layer use maths, only: pi use display, only: print_msg use utility, only: real2str ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:,:), intent(in) :: h2o_ads_reg, co2_ads_reg real(qp), intent(out) :: totmass_adsco2, totmass_adsh2o ! LOCAL VARIABLES ! --------------- integer(di) :: i, islope, l real(dp) :: thickness, slope_corr ! CODE ! ---- totmass_adsco2 = 0._qp totmass_adsh2o = 0._qp do i = 1,ngrid do islope = 1,nslope slope_corr = subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp) do l = 1,nsoil if (l == 1) then thickness = layer(l) else thickness = layer(l) - layer(l-1) end if totmass_adsco2 = totmass_adsco2 + co2_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i) totmass_adsh2o = totmass_adsh2o + h2o_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i) end do end do end do call print_msg("Total mass of CO2 adsorbed in the regolith [kg] = "//real2str(totmass_adsco2)) call print_msg("Total mass of H2O adsorbed in the regolith [kg] = "//real2str(totmass_adsh2o)) END SUBROUTINE compute_totmass_adsorbed !======================================================================= END MODULE sorption