Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ice shelf Coulomb friction law #470

Merged
merged 4 commits into from
Oct 13, 2023
Merged

Ice shelf Coulomb friction law #470

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
228976b
Added ice shelf Coulomb friction law (Schoof 2005, Gagliardini et al …
alex-huth Aug 23, 2023
9e11212
comments to describe units for basal friction parameters
alex-huth Oct 12, 2023
08d0d79
fix line length error
alex-huth Oct 12, 2023
ace1998
Merge branch 'dev/gfdl' into IS_coulomb
marshallward Oct 13, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 57 additions & 15 deletions src/ice_shelf/MOM_ice_shelf_dynamics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,11 +98,12 @@ module MOM_ice_shelf_dynamics
!! the same as G%bathyT+Z_ref, when below sea-level.
!! Sign convention: positive below sea-level, negative above.

real, pointer, dimension(:,:) :: basal_traction => NULL() !< The area integrated nonlinear part of "linearized"
!! basal stress (Pa) [R L2 T-2 ~> Pa].
real, pointer, dimension(:,:) :: basal_traction => NULL() !< The area-integrated taub_beta field
!! (m2 Pa s m-1, or kg s-1) related to the nonlinear part
!! of "linearized" basal stress (Pa) [R L3 T-1 ~> kg s-1]
!! The exact form depends on basal law exponent and/or whether flow is "hybridized" a la Goldberg 2011
real, pointer, dimension(:,:) :: C_basal_friction => NULL()!< Coefficient in sliding law tau_b = C u^(n_basal_fric),
!! units= Pa (m yr-1)-(n_basal_fric)
!! units= Pa (m s-1)^(n_basal_fric)
real, pointer, dimension(:,:) :: OD_rt => NULL() !< A running total for calculating OD_av [Z ~> m].
real, pointer, dimension(:,:) :: ground_frac_rt => NULL() !< A running total for calculating ground_frac.
real, pointer, dimension(:,:) :: OD_av => NULL() !< The time average open ocean depth [Z ~> m].
Expand Down Expand Up @@ -144,6 +145,10 @@ module MOM_ice_shelf_dynamics
real :: n_glen !< Nonlinearity exponent in Glen's Law [nondim]
real :: eps_glen_min !< Min. strain rate to avoid infinite Glen's law viscosity, [T-1 ~> s-1].
real :: n_basal_fric !< Exponent in sliding law tau_b = C u^(m_slide) [nondim]
logical :: CoulombFriction !< Use Coulomb friction law (Schoof 2005, Gagliardini et al 2007)
real :: CF_MinN !< Minimum Coulomb friction effective pressure [R L2 T-2 ~> Pa]
real :: CF_PostPeak !< Coulomb friction post peak exponent [nondim]
real :: CF_Max !< Coulomb friction maximum coefficient [nondim]
real :: density_ocean_avg !< A typical ocean density [R ~> kg m-3]. This does not affect ocean
!! circulation or thermodynamics. It is used to estimate the
!! gravitational driving force at the shelf front (until we think of
Expand Down Expand Up @@ -277,7 +282,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS)
allocate(CS%ice_visc(isd:ied,jsd:jed), source=0.0)
allocate(CS%Ee(isd:ied,jsd:jed,4), source=0.0)
allocate(CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25) ! [Pa-3 s-1]
allocate(CS%basal_traction(isd:ied,jsd:jed), source=0.0) ! [R L2 T-2 ~> Pa]
allocate(CS%basal_traction(isd:ied,jsd:jed), source=0.0) ! [R L3 T-1 ~> kg s-1]
allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10) ! [Pa (m-1 s)^n_sliding]
allocate(CS%OD_av(isd:ied,jsd:jed), source=0.0)
allocate(CS%ground_frac(isd:ied,jsd:jed), source=0.0)
Expand Down Expand Up @@ -423,6 +428,19 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
call get_param(param_file, mdl, "BASAL_FRICTION_EXP", CS%n_basal_fric, &
"Exponent in sliding law \tau_b = C u^(n_basal_fric)", &
units="none", fail_if_missing=.true.)
call get_param(param_file, mdl, "USE_COULOMB_FRICTION", CS%CoulombFriction, &
"Use Coulomb Friction Law", &
units="none", default=.false., fail_if_missing=.false.)
call get_param(param_file, mdl, "CF_MinN", CS%CF_MinN, &
"Minimum Coulomb friction effective pressure", &
units="Pa", default=1.0, scale=US%Pa_to_RL2_T2, fail_if_missing=.false.)
call get_param(param_file, mdl, "CF_PostPeak", CS%CF_PostPeak, &
"Coulomb friction post peak exponent", &
units="none", default=1.0, fail_if_missing=.false.)
call get_param(param_file, mdl, "CF_Max", CS%CF_Max, &
"Coulomb friction maximum coefficient", &
units="none", default=0.5, fail_if_missing=.false.)

call get_param(param_file, mdl, "DENSITY_ICE", CS%density_ice, &
"A typical density of ice.", units="kg m-3", default=917.0, scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "CONJUGATE_GRADIENT_TOLERANCE", CS%cg_tolerance, &
Expand Down Expand Up @@ -624,7 +642,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, &
'vi-viscosity', 'Pa m s', conversion=US%RL2_T2_to_Pa*US%Z_to_m*US%T_to_s) !vertically integrated viscosity
CS%id_taub = register_diag_field('ice_shelf_model','taub_beta',CS%diag%axesT1, Time, &
'taub', 'MPa', conversion=1e-6*US%RL2_T2_to_Pa)
'taub', 'MPa s m-1', conversion=1e-6*US%RL2_T2_to_Pa/(365.0*86400.0*US%L_T_to_m_s))
CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, &
'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m)
endif
Expand Down Expand Up @@ -720,7 +738,8 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled
real, dimension(SZDIB_(G),SZDJB_(G)) :: taud_x, taud_y !<area-averaged driving stress [R L2 T-2 ~> Pa]
real, dimension(SZDI_(G),SZDJ_(G)) :: ice_visc !< area-averaged vertically integrated ice viscosity
!! [R L2 Z T-1 ~> Pa s m]
real, dimension(SZDI_(G),SZDJ_(G)) :: basal_tr !< area-averaged basal traction [R L2 T-2 ~> Pa]
real, dimension(SZDI_(G),SZDJ_(G)) :: basal_tr !< area-averaged taub_beta field related to basal traction,
!! [R L1 T-1 ~> Pa s m-1]
integer :: iters
logical :: update_ice_vel, coupled_GL

Expand Down Expand Up @@ -2198,8 +2217,8 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask,
intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points
!! relative to sea-level [Z ~> m].
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(in) :: basal_trac !< A field related to the nonlinear part of the
!! "linearized" basal stress [R Z T-1 ~> kg m-2 s-1].
intent(in) :: basal_trac !< Area-integrated taub_beta field related to the nonlinear
!! part of the "linearized" basal stress [R L3 T-1 ~> kg s-1].

real, intent(in) :: dens_ratio !< The density of ice divided by the density
!! of seawater, nondimensional
Expand Down Expand Up @@ -2373,8 +2392,8 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac,
!! flow law [R L4 Z T-1 ~> kg m2 s-1]. The exact form
!! and units depend on the basal law exponent.
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(in) :: basal_trac !< A field related to the nonlinear part of the
!! "linearized" basal stress [R Z T-1 ~> kg m-2 s-1].
intent(in) :: basal_trac !< Area-integrated taub_beta field related to the nonlinear
!! part of the "linearized" basal stress [R L3 T-1 ~> kg s-1].

real, dimension(SZDI_(G),SZDJ_(G)), &
intent(in) :: hmask !< A mask indicating which tracer points are
Expand Down Expand Up @@ -2533,8 +2552,8 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc,
!! flow law. The exact form and units depend on the
!! basal law exponent. [R L4 Z T-1 ~> kg m2 s-1].
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(in) :: basal_trac !< A field related to the nonlinear part of the
!! "linearized" basal stress [R Z T-1 ~> kg m-2 s-1].
intent(in) :: basal_trac !< Area-integrated taub_beta field related to the nonlinear
!! part of the "linearized" basal stress [R L3 T-1 ~> kg s-1].

real, dimension(SZDI_(G),SZDJ_(G)), &
intent(in) :: float_cond !< An array indicating where the ice
Expand Down Expand Up @@ -2814,6 +2833,10 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf)
integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq
integer :: giec, gjec, gisc, gjsc, isc, jsc, iec, jec, is, js
real :: umid, vmid, unorm, eps_min ! Velocities [L T-1 ~> m s-1]
real :: alpha !Coulomb coefficient [nondim]
real :: Hf !"floatation thickness" for Coulomb friction [Z ~> m]
real :: fN !Effective pressure (ice pressure - ocean pressure) for Coulomb friction [R L2 T-2 ~> Pa]
real :: fB !for Coulomb Friction [(L T-1)^CS%CF_PostPeak ~> (m s-1)^CS%CF_PostPeak]

isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec
iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB
Expand All @@ -2825,15 +2848,34 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf)

eps_min = CS%eps_glen_min

if (CS%CoulombFriction) then
if (CS%CF_PostPeak.ne.1.0) THEN
alpha = (CS%CF_PostPeak-1.0)**(CS%CF_PostPeak-1.0) / CS%CF_PostPeak**CS%CF_PostPeak ![nondim]
else
alpha = 1.0
endif
endif

do j=jsd+1,jed
do i=isd+1,ied
if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then
umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25
vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25
unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2))
! CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction * (US%L_T_to_m_s*unorm)**(CS%n_basal_fric-1)
CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction(i,j) * (US%L_T_to_m_s*unorm)**(CS%n_basal_fric-1)
unorm = US%L_T_to_m_s*sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2))

!Coulomb friction (Schoof 2005, Gagliardini et al 2007)
if (CS%CoulombFriction) then
!Effective pressure
Hf = max(CS%density_ocean_avg * CS%bed_elev(i,j)/CS%density_ice, 0.0)
fN = max(CS%density_ice * CS%g_Earth * (ISS%h_shelf(i,j) - Hf),CS%CF_MinN)

fB = alpha * (CS%C_basal_friction(i,j) / (CS%CF_Max * fN))**(CS%CF_PostPeak/CS%n_basal_fric)
CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction(i,j) * &
unorm**(CS%n_basal_fric-1.0) / (1.0 + fB * unorm**CS%CF_PostPeak)**(CS%n_basal_fric)
else
!linear (CS%n_basal_fric=1) or "Weertman"/power-law (CS%n_basal_fric .ne. 1)
CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction(i,j) * unorm**(CS%n_basal_fric-1)
endif
endif
enddo
enddo
Expand Down