Basal hydrology scheme, new basal friction and inversion options #61
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Cleaning up some comments, e.g., changing 'Glimmer' to 'CISM'
Until now, the exponent n in the Glen flow law has been set in glimmer_physcon.F90
as an integer parameter, gn = 3.
With this commit, n is renamed n_glen (a more greppable name) for use in Glissade.
It is declared in glimmer_physcon.F90 as real(dp) with default value 3.0d0.
Since n_glen is no longer a parameter, it can now be set in the config file like other
physical 'constants' (e.g., rhoi and rhoo) that are not truly constant, but can
take different values in different models or benchmarking experiments.
To avoid changing answers in old Glide code, I retained the integer parameter gn = 3 in glimmer_physcon.F90.
This parameter is now used only in the Glide dycore (e.g., glide_velo.F90).
In Glissade, I replaced gn with n_glen in several places:
(1) In subroutine glissade_interior_dissipation_sia, the dissipation factor includes n_glen.
Note: n_glen does not appear explicitly in the 1st-order dissipation, which is proportional to tau_eff^2/efvs.
(2) In glissade_velo_sia, n_glen appears in the vertical integral for the velocity.
(3) In glissade_velo_higher, flwafact = flwa**(-1./n_glen) where flwa = A.
(4) In glissade_velo_higher, the exponent p_effstr = (1.d0 - n_glen)/n_glen is used
to compute effective_viscosity for BP, DIVA, or L1L2.
(5) In glissade_velo_higher, subroutine compute_3d_velocity_l1l2 has a factor proportional to
tau**((n_glen-1.)/2.) in the vertical integral.
For (1) and (2), n_glen was previously assumed to be an integer. Switching it to real(dp)
is answer-changing at the machine roundoff level for runs using the local SIA solver
in glissade_velo_sia.F90.
For (3), (4) and (5), n_glen replaces the equivalent expression real(gn,dp).
For this reason, answers are BFB when using the BP, DIVA or L1L2 solver.
In glissade_basal_traction, n_glen appeared in the expression for beta in the Coulomb sliding option,
HO_BABC_COULOMB_FRICTION. Here, I replaced n_glen with powerlaw_m (also = 3.0d0 by default)
to be consistent with the expressions for beta in the Schoof and Tsai laws.
In glimmer_scales, Glen's n appears in the expressions for the scaling parameters vis0 and vis_scale,
Here, I defined a local integer parameter gn = 3 so that the scales are unchanged.
It is now possible for the user to specify arbitrary values of n_glen in tests such as the slab test.
Another minor change: I added some logic to the subroutine that computes L1L2 velocities.
For which_ho_efvs = 2 = HO_EFVS_NONLINEAR, the effective viscosity (efvs) is computed from
the effective strain rate using an equation from Perego et al. (2012).
But for option 0 (efvs = constant) and option 1 (efvs = multiple of flow factor),
this strain rate equation in the code does not apply. For these options, I added an
alternative that computes velocity in terms of the strain-rate-independent efvs.
This allows us to use L1L2 for problems with constant efvs (e.g., the slab test).
The code was calling subroutine init_isostasy with isostasy = 0 = ISOSTASY_NONE. This subroutine is now called only if isostasy = 1 = ISOSTASY_COMPUTE.
I modified glissade.F90 to abort cleanly with a call to glide_finalise after an advective CFL error. This is done only when the user does *not* specify adaptive subcycling. The clean abort allows the new slabStability script to keep going, launching a new run with a shorter timestep. In subroutine glissade_flow_factor of glissade_therm.F90, I removed the FLWA_INPUT option (option 3 of whichflwa). This option is redundant with option 0, FLWA_CONST_FLWA, in which the user can set default_flwa in the parameters section of the config file, and otherwise CISM uses default_flwa = 1.0e-16 Pa^-n yr^-1. We probably should rename default_flwa to constant_flwa, but leaving it for now to avoid breaking config files in test cases. Note: This option was added by Matt Hoffman in 2014. I am unaware of test cases that use this option (most have flow_law = 0 or 2), but if there are any, we will need to fix them by switching to whichflwa = 0. In subroutine glissade_therm_driver of glissade_therm.F90, I increased the threshold for column energy conservation errors from 1.0d-8 to 1.0d-7 W/m^2. For very small timesteps of ~1.0e-6 yr, the smaller threshold can be exceeded because of roundoff errors. In subroutine glissade_check_cfl of glissade_transport.F90, I modified the fatal abort for large CFL violations (advective CFL > 10). Now, CISM aborts for CFL > 10 only when adaptive_cfl_threshold > 0, i.e. transport subcycling is enabled. This prevents excessive subcycling for egregious CFL violations. If adaptive_cfl_threshold = 0. (the default), i.e. transport subcycling is not enabled, then the code exits cleanly (with a call to glide_finalise) in glissade.F90 when advective CFL > 1. I added a TODO note to set the default value of geot (the geothermal heat flux) to 0 instead of 0.05 W/m^2. With the default value, simple prognostic tests like the dome are not mass-conserving. Not changing just yet because answers will change for the dome test.
This commit modifies the run and plot scripts for the Dukowicz slab test case, as described in Section 5 of this paper: J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a more efficient computational solution. The Cryosphere, 6, 21-34, https://doi.org/10.5194/tc-6-21-2012. The test case consists of an ice slab of uniform thickness moving down an inclined plane by a combination of sliding and shearing. Analytic Stokes and first-order velocity solutions exist for all values of Glen's exponent n >= 1. The solutions for n = 1 are derived in Dukowicz (2012), and solutions for n > 1 are derived in an unpublished manuscript by Dukowicz (2013). These solutions can be compared to those simulated by CISM. The original scripts, runSlab.py and plotSlab.py, were written by Matt Hoffman with support for n = 1. They came with warnings that the test is not supported. The test is now supported, and the scripts include some new features: * The user may specify any value of n >= 1 (not necessarily an integer). The tests assume which_ho_efvs = 2 (nonlinear viscosity) with flow_law = 0 (constant). * Physics parameters are no longer hard-coded. The user can enter the ice thickness, beta, viscosity coefficient (mu_n), and slope angle (theta) on the command line. * The user can specify time parameters dt (the dynamic time step) and nt (number of steps). The previous version did not support transient runs. * The user can specify a small thickness perturbation dh, which is added to the initial uniform thickness via random sampling from a Gaussian distribution. The perturbation will grow or decay, depending on the solver stability for given dx and dt. For n = 1, the viscosity coefficient mu_1 has a default value of 1.e6 Pa yr in the relation mu = mu_1 * eps((1-n)/n), where eps is the effective strain rate. For n > 1, the user can specify a coefficient mu_n; otherwise the run script computes mu_n such that the basal and surface speeds are nearly the same as for an n = 1 case with the mu_1 = 1.e6 Pa yr and the same values of thickness, beta, and theta. (Note: There is a subtle difference between the Dukowicz and CISM definitions of the effective strain rate; the Dukowicz value is twice as large. Later, it might be helpful to make the Dukowicz convention consistent with CISM.) I modified the plotting script, plotSlab.py, to look in the config and output files for physics parameters that are no longer hardwired. I slightly modified the analytic formulas to give the correct solution for non-integer n. This script creates two plots. The first plot shows excellent agreement between higher-order CISM solutions and the analytic solution for small values of the slope angle theta. For steep slopes, the answers diverge as expected. For the second plot, the extent of the y-axis is wrong. This remains to be fixed. I also added a new script, stabilitySlab.py, to carry out stability tests as described in: Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance of depth-integrated ice-dynamics solvers, to be submitted to The Cryosphere. The idea is that for a given set of physics parameters and stress-balance approximation (DIVA, L1L2, etc.), the script launches multiple CISM runs at a range of grid resolutions. At each grid resolution, the script determines the maximum stable time step. A run is deemed stable when the standard deviation of an initial small thickness perturbation is reduced over the course of 100 time steps. A run is unstable if the standard deviation increases or if the model aborts (usually with a CFL violation). I have run the stability script for several solvers (DIVA, L1L2, SIA, SSA) for each of two physical cases: one with thick shearing ice and one with thin sliding ice. Each suite of experiments runs in a few minutes on 4 Macbook cores for solvers other than BP. As expected, DIVA and SSA are much more stable than L1L2 and SIA. This and the previous commit correspond to the CISM code and scripts used for the initial submission by Robinson et al. (2021).
Rewrote the slab README file to describe the new command line options for runSlab.py, and the new script stabilitySlab.py.
This commit introduces new algorithms for computing the 3D velocity field in the L1L2 solver. Also, I fixed a bug so that the L1L2 solver works correctly for n_glen = 1. The goal was to improve the stability of the L1L2 solver in slab tests. At fine grid resolution (< 200 m), the solver is unstable even for very small dt, instead of following the stability limits of an SIA solver (as is the case for Yelmo). The default method (method 1) for the 3D velocity computes most quantities on the staggered grid and is unchanged. There are two new methods: * Method 2 follows the Yelmo algorithm as closely as possible, with some quantities on cell edges and some on vertices. Unlike Yelmo, there is necessarily an averaging of uvel and vvel from edges to corners at the end, since CISM assumes B-grid velocities. * Method 3 moves all intermediate calculations to cell edges, with a final averaging to vertices at the end. To support the new methods, there is a new subroutine, glissade_average_to_edges, in module glissade_grid_operators. Both new methods give results similar to method 1. Neither improves stability at very high resolution. This suggests that the stability of Yelmo might result from the overall C-grid velocity scheme rather than details of the 3D velocity computation. For each method, there is now an option to exclude the membrane stresses in the tau_xz and tau_yz terms that appear in the vertical integration factor. When these stresses are excluded, the stability curve for L1L2 solver is close to that of an SIA solver, like Yelmo.
This commit includeds the following changes: * I fixed some typos in the README file for the slab case. * In runSlab.py, I added the option to use the local SIA solver (contained in module glissade_velo_sia.F90). * I added some code which can apply an initial sinusoidal perturbation of wavelength 2*dx, instead of a random Gaussian perturbation. This is useful for getting reproducible results.
This commit adds a hybrid velocity solver as described by Robinson et al. (TC, 2021). The solver first computes SIA velocities using the local SIA solver (which_ho_approx = -1) with zero basal slip, then computes SSA velocities using the Glissade higher-order SSA solver (which_ho_approx = 1), and finally adds the two solutions. The new logic is in module glissade_velo. To use the hybrid solver, set which_ho_approx = HO_APPROX_HYBRID = 5 in the config file. For the slab test, the hybrid solver has the expected stability properties, aligned with SSA at coarse resolution and SIA at fine resolution. Answers are as expected for a dome test.
This commit streamlines the 3D velocity subroutine for L1L2 and makes a small but important change in the algorithm. The change is to evaluate the membrane stresses in tau_xz and tau_yz at layer midpoints instead of lower layer boundaries, consistent with where the SIA terms are evaluated. I found that a small vertical offset can disrupt the balance between these two terms and make the L1L2 solver unstable for the slab problem at resolutions finer than ~200 m. With the fix, the L1L2 stability curve is parallel to the SIA stability curve at fine resolution, with L1L2 being slightly more stable than SIA. I also removed the alternate L1L2 discretization methods that were added in a recent commit. The alternate strategies evaluated some terms at edges instead of vertices, and did not change the results significantly. With this commit, all terms are evaluated at either cell centers or vertices. This is the code version used for the slab tests shown for CISM in Section 3 of the paper by Robinson, Goldberg & Lipscomb (2021, TC). For these tests, I temporarily changed the energy conservation criterion in glissade_therm.F90 to avoid false non-conservation alarms when running at very fine resolution. See the comments in that module. This commit is answer-changing for L1L2, in a good way.
This commit includes a number of changes to update from python2 to python3, following the examples from other test cases. Also updated some comments, e.g. with the final citation for Robinson et al. (2022). I verified that the run, plot, and stability scripts work on Cheyenne in a python3 environment.
This commit modifies the run and plot scripts for the Dukowicz slab test case, as described in Section 5 of this paper: J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a more efficient computational solution. The Cryosphere, 6, 21-34, https://doi.org/10.5194/tc-6-21-2012. The test case consists of an ice slab of uniform thickness moving down an inclined plane by a combination of sliding and shearing. Analytic Stokes and first-order velocity solutions exist for all values of Glen's exponent n >= 1. The solutions for n = 1 are derived in Dukowicz (2012), and solutions for n > 1 are derived in an unpublished manuscript by Dukowicz (2013). These solutions can be compared to those simulated by CISM. The original scripts, runSlab.py and plotSlab.py, were written by Matt Hoffman with support for n = 1. They came with warnings that the test is not supported. The test is now supported, and the scripts include some new features: * The user may specify any value of n >= 1 (not necessarily an integer). The tests assume which_ho_efvs = 2 (nonlinear viscosity) with flow_law = 0 (constant). * Physics parameters are no longer hard-coded. The user can enter the ice thickness, beta, viscosity coefficient (mu_n), and slope angle (theta) on the command line. * The user can specify time parameters dt (the dynamic time step) and nt (number of steps). The previous version did not support transient runs. * The user can specify a small thickness perturbation dh, which is added to the initial uniform thickness via random sampling from a Gaussian distribution. The perturbation will grow or decay, depending on the solver stability for given dx and dt. For n = 1, the viscosity coefficient mu_1 has a default value of 1.e6 Pa yr in the relation mu = mu_1 * eps((1-n)/n), where eps is the effective strain rate. For n > 1, the user can specify a coefficient mu_n; otherwise the run script computes mu_n such that the basal and surface speeds are nearly the same as for an n = 1 case with the mu_1 = 1.e6 Pa yr and the same values of thickness, beta, and theta. (Note: There is a subtle difference between the Dukowicz and CISM definitions of the effective strain rate; the Dukowicz value is twice as large. Later, it might be helpful to make the Dukowicz convention consistent with CISM.) I modified the plotting script, plotSlab.py, to look in the config and output files for physics parameters that are no longer hardwired. I slightly modified the analytic formulas to give the correct solution for non-integer n. This script creates two plots. The first plot shows excellent agreement between higher-order CISM solutions and the analytic solution for small values of the slope angle theta. For steep slopes, the answers diverge as expected. For the second plot, the extent of the y-axis is wrong. This remains to be fixed. I also added a new script, stabilitySlab.py, to carry out stability tests as described in: Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance of depth-integrated ice-dynamics solvers, to be submitted to The Cryosphere. The idea is that for a given set of physics parameters and stress-balance approximation (DIVA, L1L2, etc.), the script launches multiple CISM runs at a range of grid resolutions. At each grid resolution, the script determines the maximum stable time step. A run is deemed stable when the standard deviation of an initial small thickness perturbation is reduced over the course of 100 time steps. A run is unstable if the standard deviation increases or if the model aborts (usually with a CFL violation). I have run the stability script for several solvers (DIVA, L1L2, SIA, SSA) for each of two physical cases: one with thick shearing ice and one with thin sliding ice. Each suite of experiments runs in a few minutes on 4 Macbook cores for solvers other than BP. As expected, DIVA and SSA are much more stable than L1L2 and SIA. This and the previous commit correspond to the CISM code and scripts used for the initial submission by Robinson et al. (2021).
This commit introduces new algorithms for computing the 3D velocity field in the L1L2 solver. Also, I fixed a bug so that the L1L2 solver works correctly for n_glen = 1. The goal was to improve the stability of the L1L2 solver in slab tests. At fine grid resolution (< 200 m), the solver is unstable even for very small dt, instead of following the stability limits of an SIA solver (as is the case for Yelmo). The default method (method 1) for the 3D velocity computes most quantities on the staggered grid and is unchanged. There are two new methods: * Method 2 follows the Yelmo algorithm as closely as possible, with some quantities on cell edges and some on vertices. Unlike Yelmo, there is necessarily an averaging of uvel and vvel from edges to corners at the end, since CISM assumes B-grid velocities. * Method 3 moves all intermediate calculations to cell edges, with a final averaging to vertices at the end. To support the new methods, there is a new subroutine, glissade_average_to_edges, in module glissade_grid_operators. Both new methods give results similar to method 1. Neither improves stability at very high resolution. This suggests that the stability of Yelmo might result from the overall C-grid velocity scheme rather than details of the 3D velocity computation. For each method, there is now an option to exclude the membrane stresses in the tau_xz and tau_yz terms that appear in the vertical integration factor. When these stresses are excluded, the stability curve for L1L2 solver is close to that of an SIA solver, like Yelmo.
This commit includeds the following changes: * I fixed some typos in the README file for the slab case. * In runSlab.py, I added the option to use the local SIA solver (contained in module glissade_velo_sia.F90). * I added some code which can apply an initial sinusoidal perturbation of wavelength 2*dx, instead of a random Gaussian perturbation. This is useful for getting reproducible results.
This commit streamlines the 3D velocity subroutine for L1L1 and makes a small but important change in the algorithm. The change is to evaluate the membrane stresses in tau_xz and tau_yz at layer midpoints instead of lower layer boundaries, consistent with where the SIA terms are evaluated. I found that a small vertical offset can disrupt the balance between these two terms and make the L1L2 solver unstable for the slab problem at resolutions finer than ~200 m. With the fix, the L1L2 stability curve is parallel to the SIA stability curve at fine resolution, with L1L2 being slightly more stable than SIA. I also removed the alternate L1L2 discretization methods that were added in a recent commit. The alternate strategies evaluated some terms at edges instead of vertices, and did not change the results significantly. With this commit, all terms are evaluated at either cell centers or vertices. This is the code version used for the slab tests shown for CISM in Section 3 of the paper by Robinson, Goldberg & Lipscomb (2021, TC). For these tests, I temporarily changed the energy conservation criterion in glissade_therm.F90 to avoid false non-conservation alarms when running at very fine resolution. See the comments in that module. This commit is answer-changing for L1L2, in a good way.
This is the first of a series of commits toward an improved basal water scheme for Glissade.
It is a steady-state scheme, with basal water routed conservatively down a hydraulic gradient
from input cells (where bmlt > 0) to output cells. It is based on the serial code that was written
some years ago by Jesse Johnson for Glimmer.
Although this scheme is cruder than more recent, state-of-the-art schemes such as SHAKTI
(Sommers et al., 2018), which solve time-evolving equations for the hydraulic head,
the hope is that it will allow more realistic flow simply by putting basal water
in the right locations. The current local till scheme puts water only where there
is basal melting (not in downstream locations), and is not conservative.
To use the new scheme with Glissade, set which_ho_bwat = 3 = HO_BWAT_FLUX_ROUTING in the config file.
The driver is subroutine glissade_bwat_flux_routing, in module glissade_basal_water.F90.
The scheme has 4 steps:
(1) Compute the effective pressure N for each grounded grid cell.
For now, we assume N = 0, which is equivalent to assuming that local water pressure
balances the full ice overburden pressure.
(2) Compute the hydraulic head h for each grounded grid cell. This is given by
h = z_b + p_w / (rhow*g)
where z_b = bed elevation (m)
p_w = water pressure (Pa) = p_i - N
p_i = ice overburden pressure = rhoi*g*H
N = effective pressure (Pa) = part of overburden not supported by water
H = ice thickness (m)
(3) Route basal water down the gradient of hydraulic head h. This is done by
(a) filling any depressions in h(x,y)
(b) adding small increments to h so that all water has a downhill outlet
(c) sorting the grid cells in order from low to high h
(d) initializing F = bmlt*dx*dy in each cell
(e) looping through cells from high to low h, and for each cell, partitioning
the input bwatflx among one or more down-gradient cells
For step (e), there are three routing options:
(i) D8; water is routed to the down-gradient cell with the lowest h
(ii) Dinf; water is routed to the two down-gradient cells with the lowest h,
in proportion to the gradient
(iii) FD8; water is routed to all down-gradient cells in proportion to the gradient
The original Glimmer code contained only the FD8 option. I added the others since
they are less dispersive. The user can choose the method by setting a new
config parameter, ho_flux_routing_scheme (= 0 for D8, = 1 for Dinf, = 2 for FD8).
D8 is the default.
(4) Compute the steady-state water depth based on the simple expression:
F = q * dx
where F = total water flux (m^3/s), computed in step 3
q = water flux per unit width (m^2/s), a function of grad(h) and water depth b
dx = grid cell width
Following Sommers et al. (2018, Eq. 6), we assume the flux is given by
q = (b^3 * g) / [(12*nu)(1 + omega*Re)] * -grad(h)
where b = basal water depth (m)
nu = kinematic viscosity of water (m^2/s)
Re = Reynolds number (unitless)
omega = an empirical constant (unitless)
For small Re, the flow is laminar, and for large Re, the flow is turbulent.
For now, assume laminar flow (Re -> 0).
Given F from step 3 and grad(h) from step 2, it is straightforward to solve for b.
These equations are similar, but not identical, to what is assumed in the Glimmer/Glide version.
For this reason, I did not compare try to obtain the same answers as in Glide.
To make the routing scheme more versatile, I wrote two subroutines that were not
in the old Glimmer flux-routing scheme:
(1) fix_flats (step 3b above)
This subroutine uses the algorithm of Garbrecht & Mertz (1997, J. Hydrol.)
to increment the surface elevation in flat regions, ensuring that all cells
have a downhill outlet. It repeatedly calls subroutine find_flats.
See G&M (1997) for details.
(2) find_flats
This subroutine identifies all cells without a downslope gradient.
These are regions where the surface is flat, usually after filling depressions,
and the water has no downhill outlet.
I verified that fix_flats is working, first for the example problem in G&M (1997)
and then for a dome problem with a central depression in the hydraulic head.
For diagnostics, I added output fields head and bwatflx in a new basal_hydro derived type.
I moved bwat and stagbwat to this derived type, along with some parameters that are used
for the local till model. This led to minor code changes, replacing 'temper' or 'basal_physics'
with 'basal_hydro', in several modules. Other fields and parameters could be added later
to the basal_hydro type to support new basal hydrology models, such as SHAKTI.
I coded these equations and set up a simple dome test problem with a basal melting source
beneath the dome center, where h is high. I verified that water flows downhill conservatively;
i.e., the total output flux at the margin is equal to the input flux from bmlt.
I plotted F = bwatflx for the D8, Dinf and FD8 options. As expected, FD8 gives a fairly uniform,
diffuse flow, while D8 gives a sharper flow in several distinct streams.
Dinf seems not to work well for the dome geometry because many ties are broken asymmetrically.
Next steps:
* Implement a parallel scheme on multiple tasks.
This could be done simply using gathers and scatters, or scalably using halo updates and new logic.
* Try out in more realistic ice sheet problems.
This commit includes a number of changes to allow the new flux-routing scheme
to run on multiple processors. Thus requires passing the parallel derived type
to a number of subroutines and revising the flux-routing logic.
In several iterative calculations with a stopping criterion (e.g., finding depressions,
and steps 1, 2 and 3 of the fix_flats algorithm), I replaced local with global sums.
Now, the iteration ends when a global criterion is met (e.g., no more depressions).
I modified subroutine route_basal_water as follows:
- In the serial version, the code loops from high to low cells in one iteration.
For each cell, any incoming flux is routed to neighboring downslope cells.
When we get to the lowest cell, all the water flux has reached the margin.
- In the parallel version, there are multiple iterations.
In each iteration, we loop from high to low over the locally owned cells
(not halo cells) on each processor.
The first iteration includes any water flux from basal melting (bmlt > 0).
For each local cell, the flux is routed downhill. Two things can happen:
(1) The flux reaches the low-lying margin, In this case, we are done with it.
(2) The flux is routed to a downslope halo cell. In this case, the flux in the
halo cell is communicated to the neighboring processor, and then routed
downslope to the locally owned cell adjacent to the halo.
In the next iteration, halo fluxes computed on the previous iteration are routed downhill.
When all the water has reached the margin, the iteration halts.
The total water flux is accumulated on each iteration, so that when the iteration
is done, the outgoing flux reaching the global margin should be equal to the
incoming flux received from basal melting.
To support computations of global sums, I added a parallel_global_sum interface
in the parallel modules. This could be useful in other modules too.
I also added a new halo update subroutine, parallel_halo_real8_4d, to support
efficient halo updates for arrays with two dimensions other than ewn and nsn.
Answers with a serial build are the same as for an MPI built with np = 1.
I checked that depressions are filled and flats are fixed correctly on 4 processors.
I verified that answers are the same (within roundoff) on 1 versus 4 processors
for a dome problem with a simple hydraulic head field that has flow across processors.
It might be possible to make the routing algorithm more scalable, e.g. by reducing the
number of global sums. However, this might not be worth the effort, if the
flux-routing scheme remains much cheaper than the velocity solver.
This commit includes some bug fixes and other changes in the parallel flux-routing scheme:
* Fixed a bug in the Dinf routing option. Results now look sensible.
* Do not apply bmlt < 0 (i.e., refreezing) to the basal water flux.
To conserve water, refreezing will need to be handled separately.
* Changed the criteria for bwat_mask = 0. This is the mask that identifies cells
through which basal water can flow. We do not want to exclude ice-free cells
in the interior, but we want to count all water that leaves the ice sheet.
Now, bwat_mask is set to 1 in the following cells:
(a) floating or open ocean
(b) cells at the edge of the global domain. I added a subroutine,
parallel_global_edge_mask, to identify these cells at initialization.
(c) ice-free cells with overwrite_acab_mask = 1
* Added an option for overwrite_acab: OVERWRITE_ACAB_INPUT_MASK = 3.
With this option, cells where the SMB is overwritten (usually with a negative value)
are read from the input file. For the dome hydro problem, this mask can be used
to zero out basal water outside the original dome boundary.
* Increased count_max, the max number of iterations for the flux routing, to 50.
A new iteration is needed whenever there is a nonzero flux in one of more halo cells.
With Dinf or FD8, there can be multiple crossings of the boundary between processors
as water flows down the gradient. (Up to ~40 on a 4km Greenland grid.)
This commit contains several minor changes to support the flux-routing hydrology scheme in runs with Northern Hemisphere paleo ice sheets * Added a new effective pressure option: which_ho_effecpress = HO_EFFECPRESS_BWAT_RAMP = 5. This is similar to the existing option HO_EFFECPRESS_BWAT = 4. The only difference is that as bwat increases from 0 to bwat_till_max, the effective pressure ramps down linearly, unlike the more complex formulation of Bueler and van Pelt (2015). * Added a new logical option, smooth_input_usrf. When this option is set to true, the initial upper surface elevation field (usrf) is smoothed, and the thickness is then adjusted to be consistent. This is helpful for a 4-km N. Hemisphere simulation using the input file cism_USGS-huy3-S1D_4km.nc. In this file, parts of the GrIS have rough topography and fairly smooth thck, leading to rough usrf and large surface gradients that crash the solver in the first few time steps. With several smoothing passes, the code starts cleanly. The smoothing is done in subroutine glissade_smooth_usrf, in glissade_utils.F90. Note that usrf is smoothed only for grounded ice, and not for floating ice, ice-free ocean, or ice-free land. * Made the fill_depression routine more efficient by restructuring with an outer and inner loop. Halo updates are called only from the outer loop. Although more efficient than before, it is still very expensive to fill depressions with the current algorithm. The cost of the entire code more than doubles on a 4-km N. Hemisphere grid. In a future commit, I will try to implement a different algorithm that scales better. * Added the max value of bmlt_ground in the diagnostic log file
This commit adds an efficient algorithm for filling depressions in the head field
when running the flux-routing hydrology scheme.
The old scheme was very slow to converge on large grids such as the 4km
Northern Hemisphere grid.
The new scheme is based on the algorithm of Planchon and Darboux (2001).
The basic idea is:
* Initially, set phi = phi_in on the boundary, and set phi = a large number elsewhere
(where phi is the head field).
* Loop through the domain. For each cell c, with value phi(c) not yet fixed to a known value,
compute phi_min8(n), the current minimum of phi in the 8 neighbor cells.
- If phi_in(c) > phi_min8(n) + eps, then set phi(c) = phi_in(c) and mark that cell as having a known value,
since phi(c) cannot go any lower.
- If phi_in(c) < phi_min8(n) + eps, but phi(c) > phi_min8(c) + eps, set phi(c) = phi_min8(n) + eps.
Do not mark the cell as having a known value, because it might be lowered further.
* Continue until no further lowering of phi is possible. At that point, phi = phi_out.
Here, eps is a small number greater than zero, which ensures that there are no flat surfaces
when the depression-filling is done. Thus, it is no longer necessary to call fix_flats.
On the 4km N.H. grid, the number of depression-fill iterations is reduced from several hundred
per time step to ~10.
I ran the EISMINT-2 tests with Glide for the first time in recent memory
to see if they were still working. They were not. Experiment G, which
has substantial basal sliding, crashed with a checkerboard instability
in the thickness field.
The main problem was a sign error in glide_velo.F90, subroutine slipvelo.
The code originally contained this line:
model%velowk%fslip = rhograv * btrc
which I had replaced some time ago with this line:
model%velowk%fslip = rhoi * grav * btrc
However, the correct replacement includes a minus sign:
model%velowk%fslip = -rhoi * grav * btrc
There are a number of constants and fields in this module that are defined as negative.
I added some comments to clarify which expressions implicitly include minus signs.
A minor fix: I had added bmlt_ground as a diagnostic field without making sure
it is allocated for Glide runs. Now, bmlt_ground is always allocated.
For Glide, which has no melting beneath floating ice, bmlt_ground is simply
a copy of bmlt.
To each config file, I added diagnostics at the central dome point (31,31) every 5000 years.
This makes it possible to read relevant results (e.g., total area and volume, and
basal temperature at the divide) directly from the log file.
I changed bmlt_ground back to bmlt in the list of output fields.
I then ran the full suite of EISMINT-2 experiments with Glide. All runs finished cleanly.
I compared results from Experiments A, B, C, D, G and H with those in the paper
by Payne et al. (2000, J. Glaciol.).
Glide results are generally consistent with the published model means.
The paper does not say which specific model was Tony's early version of Glide,
but it probably was V, W or Z based on the spoke patterns in Fig. 11 for Expt H.
Glide has changed enough that we would not expect to duplicate his answers.
This commit is answer-changing (in a good way) for Glide runs with nonzero basal sliding.
This commit adds a directory with config files for running several EISMINT-2 experiments
(A, B, C, D, F, G and H) with the Glissade local shallow-ice solver.
These are the seven experiments discussed in the paper by Payne et al. (2000, J.Glaciol.).
Tests I, J, K and L are omitted.
The Glissade local SIA solver is similar to Glide, except that the continuity equation
is solved using an explicit upwind-weighted advection scheme (usually incremental remapping)
instead of an implicit diffusion scheme. See glissade_velo_sia.F90 for the solver algorithm.
The different settings, compared to the Glide tests, are as follows:
* dycore = 2 (Glissade)
* evolution = 3 (IR transport)
* vertical_evolution = 0 (not constrained by upper kinematic BC)
* which_ho_approx = -1 (Glissade local SIA solver)
* which_ho_sparse = 3
* ice_limit = 1 (minimum thickess for active ice = 1 m)
Note: The Glissade local SIA scheme does not need a sparse solver,
but the default value (which_ho_sparse = 0) triggers an error in parallel runs.
Tests are typically run as follows:
> ./cism_driver e2.a.config
where cism_driver is a symbolic link to the executable:
cism_driver@ -> ../../../builds/mac-gnu/cism_driver/cism_driver
The Glissade SIA solver is slower than Glide.
Unlike Glide, however, the Glissade solver can be run in parallel, e.g.
> mpirun -n 4 ./cism_driver e2.a.config
This commit modifies module eismint_forcing.F90 so that the global temperature
and mass balance forcing are computed correctly on multiple processors.
I ran experiments A and G to completion and compared to Glide results.
The results are different (e.g., Glissade gives a higher dome), but
are in the same ballpark as the results in Payne et al. (2000).
One config file, e2.a.config.diva, is included with typical higher-order settings
for the Glissade DIVA solver. Other config files could be modified in a similar way.
In some runs, a spurious nonzero tf_anomaly was appearing in subroutine glissade_bmlt_float_thermal_forcing, as a result of not being initialized to zero during each pass through the calling subroutine in glissade.F90. In the first few timesteps, the anomaly grew large enough to trigger a fatal error in the check for reasonable values of thermal forcing. This commit adds logic to initialize tf_anomaly and tf_anomaly_basin correctly. It's unclear why this bug didn't show up before.
This is an initial implementation of the sliding law suggested by Zoet & Iversion (2020, Science).
This law takes the form (see Eq. 3 in the ZI paper):
tau_b = C_c * N * [u_b/(u_b + u_t)]^(1/m), Eq. 3 in ZI(2020)
where C_c = a constant in the range [0,1], aka coulomb_C
N = effective pressure
u_t = threshold speed controlling the transition between powerlaw and Coulomb behavior
m = powerlaw exponent
Here, C_c and m correspond to tan(phi) and p, respectively, in ZI Eq. 3.
This law is implemented as which_ho_babc option HO_BABC_ZOET_IVERSON = 7,
replacing the unsupported option HO_BABC_YIELD_NEWTON.
C_c can be implemented either as a constant (default value = 0.42) or as a 2D field (coulomb_c_inversion).
For the latter, I implemented an inversion procedure, closely following the procedure for Cp inversion
with the Schoof and powerlaw options.
The threshold speed u_t is a new config parameter.
Tim van den Akker already added ZI on a branch created from main.
This commit replicates Tim's changes on a branch that is rebased onto the recent basal hydrology changes,
and adds code changes for C_c inversion.
Note: The new code has not yet been tested thoroughly.
This commit may be squashed with other commits after debugging.
This commit includes a number of changes in the options for computing the
yield stress in Coulomb-type friction laws. These are laws in which the
yield stress is written as N*C_c or N*tan(phi), where N is effective pressure,
C_c is a scalar in the range [0,1], and phi is a friction angle, such that
tan(phi) is in the range [0,1].
Note that C_c and tan(phi) are functionally equivalent, and CISM generally
uses the C_c terminology (coulomb_c in the code).
The goal is to allow mixing and matching of options for specifying or reducing
the yield stress: e.g. via ocean connection, or the presence of meltwater at the bed,
or permanently weak till.
The main changes are:
* There is a new set of options for which_ho_effecpress; some are removed and some are added.
* The ocean connection option is now applied simply by setting p_ocean_penetration > 0,
without needing to set which_ho_effecpress = 3.
* There is a new option called which_ho_coulomb_c for setting till strength.
* There is a new option which_ho_powerlaw_c, analogous to which_ho_coulomb_c.
* The options which_ho_cc_inversion and which_ho_cp_inversion are deprecated.
Instead, the user sets which_ho_coulomb_c/powerlaw_c = 1 to compute these
coefficients by inversion, and which_ho_coulomb_c/powerlaw_c = 2 to read these
coefficients from an external file (which may have been obtained by inversion).
* Added 2D fields coulomb_c_2d and powerlaw_c_2d to go with these new options,
superseding coulomb_c_inversion and powerlaw_c_inversion.
The new fields are part of the basal_physics derived type (instead of the inversion type).
With these changes, config files used for ISMIP6 experiments will not, in general,
run cleanly and give the same result with the new options. Several changes will be needed
in the CESM namelists.
More details are given below.
*****
The new options for which_ho_effecpress are:
[0] N is equal to overburden pressure, rhoi * g * H.
[1] N is reduced where the bed is at or near pressure melting point, with a linear ramp.
[2] N is reduced where basal water is present (bwat > 0), with a linear ramp.
[3] N is reduced where there is meltwater flux at the bed (bwatflx > 0), with a linear ramp.
[4] N is reduced where basal water is present (bwat > 0), following Bueler & van Pelt.
Options 0, 1, and 4 are the same as current options.
Option 2 is the same as the old option 5.
Option 3 is new, replacing the old ocean_penetration option.
The options for which_ho_coulomb_c are:
[0] C_c is equal to a spatially uniform constant, the parameter coulomb_c.
I changed the default value from 0.42 to 1.0. In other words, yield stress
is equal to overburden pressure by default, but can be overridden by
setting coulomb_c = 0.50, or a similar value, in the config file.
[1] C_c is found by inversion, nudging the simulated ice thickness toward a target thickness.
[2] C_c is a 2D field read from an external file. The values could be based on geology
or obtained from a previous inversion.
[3] C_c is a function of bed elevation, following the pseudo-plastic law as typically applied.
Note: Option 3 is not supported with this commit; will add support in an upcoming commit.
The options for which_ho_powerlaw_c are:
[0] C_p is equal to a spatially uniform constant, the parameter powerlaw_c.
The default value remains 1.0d4 Pa m^(-1/3) yr^(1/3), which is appropriate
for a power law with exponent m = 3.
[1] C_p is found by inversion, nudging the simulated ice thickness toward a target thickness.
[2] C_p is a 2D field read from an external data set. The values could be based on geology
or obtained from a previous inversion.
There is no analog to which_ho_coulomb_c = 3.
A conceptual point: The basal melt mechanisms that decrease N by increasing water pressure
could equally be viewed as mechanisms that decrease C_c by weakening the till.
Here, we view these mechanisms as decreasing N. Thus, N becomes a dynamic field
that can vary rapidly with changes in the basal water system or ice thickness, whereas
C_c is viewed as a slowly varying field embodying quasi-permanent properties of the till.
With which_ho_coulomb_c = 2 and isostasy turned on, C_c will vary with the bed elevation,
but this is generally a slow change.
If the user mistakenly sets which_ho_effecpress = 3 to specify ocean penetration with p > 0,
the code should get the same answer as before, provided bwatflx = 0.
But it is safer to set which_ho_effecpress = 0 in this case, while setting p > 0.
Another small change: Changed the units of bwatflx from m^3/yr to m/yr, so that this flux
is not proportional to the grid cell area.
I set up test cases for ISMIP6-style Antarctic spin-up runs and confirmed that the answers
are BFB with the appropriate choices of new options.
This commit implements option which_ho_coulomb_c = HO_COULOMB_C_ELEVATION = 3, with coulomb_c varying linearly between a max value at high bed elevations (>= bedmax) to a min value at low bed elevations (<= bedmin). This option generalizes the procedure used in the pseudo-plastic basal friction law to set tan(phi), which is equivalent to coulomb_c. For details, see the git log for the previous commit. To use the new option, users should set the config parameters coulomb_c_max, coulomb_c_min, coulomb_c_bedmax, and coulomb_c_bedmin, if not using the default values. The new coulomb_c option can be applied to any basal friction law with a coulomb_c term. Currently, these are friction laws 2 (new pseudo-plastic), 7 (Zoet-Iverson), 10 (basic Coulomb), 11 (Schoof), and 12 (Tsai). Previously, bed elevation dependence was hardwired for the pseudo-plastic friction law (the default in CESM Greenland runs), but was not supported for other friction laws. For backward compatiblity, CISM now supports two flavors of pseudo-plastic law: * HO_BABC_PSEUDO_PLASTIC_OLD = 3, the old option, with hardwired elevation dependence * HO_BABC_PSEUDO_PLASTIC = 2, the new option, with optional elevation dependence The older option can be removed when it is no longer the CESM default. To make room for the new option as #2, I set the little-used HO_BABC_YIELD_PICARD = 15. To streamline the code, powerlaw_c_2d and coulomb_c_2d are now computed near the start of subroutine calcbeta, instead of being computed separately for each friction law. Note: While adding the new option, I found that the old option has a minor bug: tanphi is a function of topg(ew,ns). However, tanphi should be colocated with beta on the staggered grid, whereas topg lives on the unstaggered grid. I left the bug in place for backward compatibility, given that this option will be deprecated going forward. Note: The default value of coulomb_c_min is 1.0e-3. For the pseudo-plastic law, a more appropriate value is coulomb_c_min ~ 0.10, which is close to tan(phimin) with phimin = 5 degrees. In runs with a pseudoplastic law, answers are BFB with which_ho_babc = 3, and are modestly different (as expected) with which_ho_babc = 2 combined with which_ho_coulomb_c = 3. In Antarctic runs with other friction laws, answers are BFB except when using the new options.
Previously, basal_physics%powerlaw_c and basal_physics%coulomb_c were constants, and basal_physics%powerlaw_c_2d and basal_physics%coulomb_c_2d were 2D fields. The usual CISM convention is not to put '2d' in the names of fields except to distinguish them from 3d fields (as in the case of uvel_2d, vvel_2d). Instead, we put 'const' or 'constant' in the names of parameters that need to be distinguished from 2D fields. With this commit, the constants are named basal_physics%powerlaw_c_const and basal_physics%coulomb_c_const. The 2D fields are simply basal_physics%powerlaw_c and basal_physics%coulomb_c. This commit is BFB. However, it requires changing config files. In the parameters section, powerlaw_c -> powerlaw_c_const and coulomb_c -> coulomb_c_const. In lists of output fields, powerlaw_c_2d (or the older powerlaw_c_inversion) becomes powerlaw_c, and similarly for coulomb_c.
When running ISMIP6 experiments, it was convenient to read in some parameters and fields from the input file. However, this has led to some confusion and errors. In particular, when gamma0 can either be read from an input or restart file or specified in the config file, there is a risk of overwriting the desired value. With this commit, gamma0 is no longer an I/O option. It is not listed in glide_vars.def or written to restart files. The only way to set gamma0 /= 0.0 (the default is 0.0) is to specify a nonzero gamma0 in the config file. Thus, gamma0 is now treated as most other parameters are treated. I also removed the bmlt_float_ismip6_magnitude option. This option was used to specify whether CISM was running with low-, median-, or high-sensitivity gamma0, as defined for ISMIP6, and thereby choose specific values of gamma0 and deltaT_basin. Six deltaT_basin variants (deltaT_basin_nonlocal_median, deltaT_basin_local_pct5, etc.) are no longer in the code. If present, deltaT_basin is still read from the input file, but the user must put the correct field in the input file, instead of letting CISM choose from among six possible fields. When inverting for deltaT_basin, there is no need to put it in the input file; it is initialized to zero everywhere. Finally, I removed the thermal_forcing_baseline field, which was no longer used but was simply copied to thermal_forcing. This commit is BFB, apart from removing some old options. Note: Some older Antarctic input files contain a field called thermal_forcing_baseline. This should be changed to thermal_forcing, else it will not be read in.
This commit removes most of the plume model that was developed several years ago but was never robust. The plume model was written to support option whichbmlt_float = BMLT_FLOAT_MISOMIP = 5, with the goal of computing sub-shelf melt rates given the prescribed MISOMIP temperature and salinity profiles. I left the option in the code, since it might be supported with a plume or reduced-order ocean model that is still to be developed. However, the code now aborts with a fatal error if the user selects this option. In glissade_bmlt_float.F90, I removed most of the plume code but kept a few subroutines that might be useful later.
This commit removes the old bmlt inversion option, which_ho_bmlt_inversion, which inverts for bmlt_float in floating grid cells, adjusting the melt rate to bring the thickness in each grid cell closer to an observational target. This scheme was used for the original ISMIP6 Antarctic runs but was never very robust. The melt rate is sensitive to the ice thickness, so it is hard to converge on stable values. The resulting melt rate field is noisier than real melt rates would be. The bmlt_basin inversion option (which_ho_bmlt_basin_inversion) remains. This scheme inverts for deltaT_basin in each basin and is more robust.
The steady-state flux-routing scheme computes the basal water flux (bwatflx) based on the assumption that the total water inflow from basal melting is equal to the outflow to the ocean. It then diagnoses bwat from bwatflx, given an expression for the water velocity in terms of grad(head). Previously, the effective pressure could be reduced in subroutine calc_effecpress based on either bwat or bwatflx. However, the bwat diagnosed in the flux-routing scheme is problematic for basal thermodynamics. The thermal solver assumes that the bed is at the melting point wherever bwat > 0. As a result, a very small bwat beneath frozen ice in the flux router can drive large, abrupt increases in basal temperature. The solution is to diagnose bwat locally in the flux router, without returning bwat to glissade. As a result, the option which_ho_effecpress = HO_EFFECPRESS_BWAT = 2 can no longer be used with the flux router. It should be used only with the local till model, which prognoses bwat. To reduce N when using the flux router, set which_ho_effecpress = HO_EFFECPRESS_BWATFLX = 3. For consistency, effecpress_bwat_threshold is now used for option 4 (Bueler-Van Pelt) as well as option 2 (linear ramp). The BvP scheme has a related parameter, bwat_till_max, which normally should be set to the same value as effecpress_bwat_threshold. Both parameters now have defaults of 2 m. I verified that a sample Greenland run with the pseudo-plastic law and local till is BFB. A minor change: I commented out the call to subroutine effective_pressure in the flux router, and added comments explaining that we assume N = 0 for the purpose of computing the head. With these changes, Antarctic simulations are more stable than before. However, with which_ho_effecpress = 3 and effecpress_bwat_threshold > 1 m/yr, only the largest Antarctic channels have significant reductions in N. We should consider replacing the linear ramp with other functional forms.
The flux-routing basal water scheme returns a 2D bwatflx field, which is used to reduce the effective pressure N. Previously, we used a linear ramp function to compute N from bwatflx. However, bwatflx can range over several orders of magnitude, and the marginal reduction in N per unit of added bwatflx could decrease gradually as bwatflx grows. With this commit, N is prognosed as a function of Fw/F0, where Fw = bwatflx and F0 is an empirical constant. More precisely, the prognosed quantity is called f_effecpress (or f_N for short), defined as the ratio of effective pressure to overburden pressure, N/N_o = N/(rhoi*g*H). That is, f_N is the fraction of overburden that is *not* balanced by water pressure. It lies in the range (0,1] and evolves as follows: df_N/dt = [1 - f_N*(Fw/F0)] / tau_N, where tau_N is a relaxation timescale. With 'f_N' in the numerator of the rhs, it becomes harder to reduce N as f_N approaches zero. The '1' in the numerator nudges f_N back toward 1 when Fw is small or zero. The asymptotic value attained with a given Fw is f_N = min(F0/Fw, 1). Note: Initially, I tried diagnosing N based on the instantaneous value of (Fw/F0), but this can lead to unstable oscillations in ice speed and thickness. For now, the default values of the new constants are tau_N = 100 yr and F0 = 0.01 m/yr. These can be set in the [parameters] section of the config file. I made a related change in the ZI law. This law has the form tau_b = N * C_c * [u_b/(u_b + u_t)]^(1/m). I replaced 'N' with min(N, N_max), where N is the actual effective pressure and N_max is another empirical constant. The idea is that for thick, slow-sliding ice in the power-law regime, the sliding coefficient (i.e., the term that multiplies u_b^(1/m)) can asymptote instead of continuing to increase with higher N. This makes the ZI law more like the Schoof law in the power-law limit. The default N_max = 1.e8 Pa, an overburden corresponding to an unrealistic ~10 km of ice. Thus, N is not capped in the ZI law by default, but only when the user sets a smaller N_max in the config file. Finally, I changed the way N is computed when it can be reduced by either basal melting (which_ho_effecpress >= 1) or by p > 0. Previously, N could be multiplied successively by two numbers < 1 when both effects were active. On further reflection, one reduction should not necessarily reinforce the other. The new computation is N = min(N(bwatflx), N(ocean_p)). That is, we choose the minimum value from the two calculations.
This commit fixes a bug in the computation of effective pressure from bwatflx. A halo update of bwatflx was missing, leading to an error in staggered effecpress and resulting ice speeds at processor boundaries. Thanks to Sarah Bradley for spotting the error. The fix is to add a halo update for bwatflx at the end of the flux-routing scheme, and/or a halo update for effecpress before computing effecpress_stag. Either update fixes the bug, but to be on the safe side, I added both updates. I also added and revised some basal water and effecpress diagnostic prints.
I added logic to prevent evaluating a term with thck_gradient_ramp in the denominator unless thck_gradient_ramp > 0. This turned up in recent debugging; not sure why it didn't turn up earlier.
Until now, the inversion for powerlaw_c and coulomb_c has used an
observational thickness target (or equivalently, a surface elevation target).
With this commit, it is possible to invert for powerlaw_c or coulomb_c based on
a thickness target, a velocity target (observed surface speed), or both.
The two targets are combined in one subroutine, with similar logic for each target.
The rate of change of C_p or C_c is proportional to:
(H - Hobs) / babc_thck_scale with a thickness target, and/or
(v - vobs) / babc_velo_scale with a velocity target.
The rate is damped by a term proportional to dH/dt for a thickness target,
but there is no damping proportional to dv/dt. (I tried adding a damping term
proportional to dv/dt, but it led to oscillations in coulomb_c.)
To support the velocity target, these parameters and fields have been added:
babc_velo_scale = scale for velocity differences
usfc_obs = x component of observed surface speed
vsfc_obs = y component of observed surface speed
velo_sfc_obs = observed surface speed = sqrt(usfc_obs^2 + vsfc_obs^2)
Note: babc_thck_scale = babc_velo_scale = 0.0 by default.
Previously (e.g., for ISMIP6 runs), babc_thck_scale was equal to 100 m by default.
To enable inversion based on thickness, velocity, or both at once,
the user can set either or both parameters to a positive value in the config file.
For runs with velocity inversion, I am using a value of 200 m/yr.
To reduce code duplication, I consolidated several subroutines in glissade_inversion.F90.
The same subroutines now handle both powerlaw_c and coulomb_c inversion,
with which_ho_powerlaw_c and which_ho_coulomb_c determining the appropriate arguments
to pass between subroutines.
I verified that answers are BFB with the consolidation.
I removed dthck_dt from the list of inversion restart variables.
It is not, in fact, needed for exact restart.
The scheme seems to be working as desired in a 3 ka Antarctic spin-up.
Still need to do a close analysis of runs with and without a velocity target.
One challenge in Antarctic spin-ups is to prevent the retreat of the East Thwaites grounding line early in the run. The ice can thin by several hundred meters and never recover. With p = 1, the entire Thwaites basin can collapse. With this commit, users can phase in the effect of ocean_p (aka p_ocean_penetration or simply 'p') on effective pressure N by specifying a relaxation timescale > 0. The timescale is a new parameter, basal_physics%ocean_p_timescale. When the timescale > 0, the effective pressure is relaxed using a new prognostic field, basal_physics%f_effecpress_ocean_p. This is a scalar field in the range [0,1] that specifies the fractional value of N relative to overburden. If p > 0, then N is multiplied by f_effecpress_ocean_p for each gridcell. Typically, the original N is the overburden pressure rhoi*g*H, unless N is also reduced by meltwater at the bed (e.g., with a hydrology model). With a timescale of zero (the default value), N = rhoi*g*H*(1 - Hf/H)^p. In this case, f_effecpress_ocean_p = (1 - Hf/H)^p (where Hf = flotation thickness). With a timescale > 0, we take (1 - Hf/H)^p as a target fraction, but relax toward that value over a period of ocean_p_timescale. As a result, N does not immediately drop to a low value in regions where Hf is close to H. This gives the ice time to thicken, so that Hf/H becomes smaller and the target fraction is larger. Thus, f_effecpress_ocean_p may never reach the low value that was the initial target. There was already a field called f_effecpress that acts in a similar way by phasing in the contribution of bwatflx to reducing N. This field is now call f_effecpress_bwat to avoid confusion with f_effecpress_ocean_p. In new Antarctic simulations, I set basal_physics%ocean_p_timescale = 500 yr, the same value as inversion%babc_timescale. This change inhibits retreat of the E. Thwaites GL (there is still some thinning, but much less than before) and generally reduces thickness biases for lightly grounded ice. In making this change, I reverted to an earlier method for compounding reductions in N. This method is multiplicative. Thus, if both basal water and ocean_p act to reduce N, they act in concert; we can multiply N by f_effecpress_bwat and again by f_effecpress_ocean_p. The more recent method was to choose the minimum of (1) N as reduced by basal water and (2) N as reduced by ocean_p > 0. Going back to the earlier method allows lower N and reduces the need for low C_c in some regions. Note: Setting ocean_p_timescale = 0 (the default) results in roundoff-level answer changes compared to previous code, because some operations have been reordered. Other changes: - Updated the ho_whicheffecpress descriptions in glide_setup. - Modified the interface for glide_define_restart_variables to pass in 'model' instead of 'model%options'. This allows restart variables to be added based on whether certain parameters (e.g., p_ocean_penetration) are zero or nonzero. - Changed a comment in glissade_bmlt_float
When inverting for bmlt_basin, CISM now sets the area and volume targets based on all the floating ice in the basin, plus the lightly grounded ice for which f_flotation > -(bmlt_basin_flotation_threshold). In the previous logic, floating ice was included in the target only where f_flotation < bmlt_basin_flotation_threshold. The goal is to more accurately simulate ice-shelf thickness by including relatively more floating ice compared to grounded ice in the target. In the latest runs, I set bmlt_basin_flotation_threshold = 200 m instead of 500 m, to exclude ice that is fairly well grounded. Also, when initializing coulomb_c inversion, coulomb_c is set to coulomb_c_const instead of 1.0. Assuming coulomb_c_const ~0.5, this allows ice to flow a bit faster at the start of the run without leading to instability.
Until now, the mask-based calving scheme has assumed that for all floating ice cells with calving_mask = 1, any ice in the cell is removed immediately. The calving front cannot advance, even if the flux at the CF is very large. This can make it harder to spin up ice shelves accurately in regions where the ice is (or recently was) grounded near the CF, as for E. Thwaites. With this commit, it is possible to remove ice in masked cells incrementally, by specifying a positive calving timescale in the config file. For instance, if calving_timescale (aka calving%timescale) is 10 years, then we remove 1/10 of the thickness each year in cells with calving_mask = 1. This, ice that starts with a thickness of 300 m would have a thickness of 270 m after 1 year, and roughly 1/e of its initial thickness after 10 years. To prevent thin ice from lingering indefinitely, this thinning should be supplemented with thickness-based calving, so that ice thinner than a given threshold, calving%minthck, is removed rapidly. To run with this new option, the user sets marine_margin = 5 (the usual setting for mask-based calving) in the config file, along with positive values of calving_timescale and calving_minthck. A large timescale allows the calving front to advance well beyond the observed front. The goal is to choose a timescale that allows the CF to advance modestly, possibly with grounding on ridges, where favored by the dynamics. I have done some Antarctic runs with calving_timescale = 10 yr and calving_minthck = 100 m. This seems to have a beneficial effect on E. Thwaites, but more testing is needed to find optimal parameters. I also added logic to set calving_mask = 0 (i.e., do not calve floating ice) not only where thck > 0 in the input file, but also where usfc_obs or vsfc_obs > 0 in the input file. Previously, calving_mask was based on the input thickness alone. The goal of this change is to let ice shelves (e.g., Thwaites) advance to grid cells where ice has existed during the observational period, even if thck = 0 in the latest data. This includes parts of the Thwaites shelf that are ice-free in BedMachine data, but have nonzero speeds in earlier velocity data.
In recent Antarctic spin-ups with inversion for deltaT_basin, I was unable
to get a steady-state grounding line near the observed position in the central
Ronne ice shelf, aka the Southern Weddell Sea Embayment. The GL there is perched
on the Bungenstock Ice Rise (BIR) and is easily dislodged.
In spin-ups, the GL first advances, leading to an increase in deltaT_basin
by ~0.2 degrees. This results in excessive melting, driving back the GL.
Instead of settling on the BIR, the GL recedes past the Robin Subglacial Basin,
leading to significant SLR compared to the present day.
We expect that this basin might collapse under future scenarios,
but we don't want it to collapse during the spin-up.
The alternative implemented here is to replace the deltaT_basin inversion
with inversion for a basin-scale flow factor. Until now, there has been one factor
(called 'flow_factor' in the config file) for grounded ice, and another factor
('flow_factor_float) for floating ice. Typically, we set flow_factor_float to ~0.4
so the large ice shelves will not flow faster than observed.
The goal with the new inversion is to vary flow_factor_float in each basin.
Where the basin average is too thick, flow_factor_float is increased,
making the ice softer. Where the basin average thickness is too thin,
flow_factor_float is decreased, making the ice stiffer. In this way we hope
to match observed shelf thickness without imposing unrealistic melt rates.
The main code changes are as follows:
* New 2D I/O field, flow_factor_basin, which is initialized to a constant and then
adjusted for each basin using a prognostic equation similar to that for deltaT_basin
* New config option, which_ho_flow_factor_basin, with three choices: (0) set to a constant,
(1) obtain by inversion, (2) read from an external file (based on previous inversion).
* New config parameters, flow_factor_basin_thck_scale and flow_factor_basin_timescale,
used in the prognostic equation for flow_factor_basin. Defaults are 100 m and 500 yr.
* New subroutine, glissade_inversion_flow_factor_basin, in the glissade_inversion module.
Like the deltaT_basin subroutine, the new routine is called during the diagnostic solve.
This subroutine prognoses flow_factor_basin.
For now, this factor is forced to stay in the range [0.20, 3.0].
* Subroutine glissade_flow_factor, which computes flwa (aka 'A') in the flow law,
now has a 2D field, flow_enhancement_factor_float, as an input argument.
This replaces the constant parameter flow_factor_float.
* Moved the basin_sum and basin_average subroutines to module glissade_utils
* Changed the initialization of coulomb_c. Now, coulomb_c is initialized to
coulomb_c_const for grounded cells and coulomb_c_min for floating and ocean cells.
Thus, when the grounding line advances, it encounters low basal friction and
is less likely to continue advancing.
* Option which_ho_bmlt_basin_inversion is now called which_ho_bmlt_basin, with four choices:
(0) set deltaT_basin = 0 for all basins, (1) obtain by inversion, (2) read from file,
and (3) apply values ISMIP6 values obtained by inversion.
Option (3) hardcodes the values obtained by N. Jourdain et al. (2019), supplemented by
values Nico computed for Lipscomb et al. (2021). I added Nico's deltaT_basin values
in subroutine glissade_bmlt_float_thermal_forcing_init to support this option.
* Config parameter inversion_bmlt_basin_flotation threshold is now called
inversion_basin_flotation_threshold (without 'bmlt'), and similarly for
inversion_basin_mass_correction and inversion_basin_number_mass_correction.
These parameters apply to either basin inversion option.
I compared options (0) and (3) in spin-ups. In some basins, one or the other looks better,
but neither seems superior overall. For example, Nico's negative adjustment in basin 11
prevents the George VI shelf from thinning too much, whereas the positive adjustment
in basin 9 leads to excessive Amundsen melting.
This suggests we might as well set deltaT_basin = 0 everywhere as the default,
using viscosity to correct for thickness biases.
I ran a number of spin-ups to test the new flow_factor_basin inversion scheme.
The scheme is working as desired. It turns out that the majority of basins
reach either their max value (3.0) or min value (0.20).
However, the Ross and Ronne basins typically get a value between 0.5 and 2.0.
A relatively small change in deltaT_basin can drive relatively large changes
in basal melt rates and can thereby lead to a large change in flow_factor_basin,
especially for smaller shelves.
The hope is that flow_factor inversion will make the grounding line more stable
in the BIR region of the Ronne shelf. However, the BIR GL still retreats under many
parameter choices, including p = 0.25 or greater with the standard ISMIP6
16-basin scheme. I tested a 28-basin scheme that isolates the BIR/Institute/Moeller region
as a single basin. We can then obtain a flow factor lower than in the neighboring
Filchner-Ronne basins, which tend to be thicker than observed. We will continue these tests.
An unrelated fix: Subroutine calc_effecpress is no longer called during the first
diagnostic solve after a restart, since this can break exact restart.
The prognostic equation when inverting for either C_p or C_c has been:
(1/C) * dC/dt = -(1/tau) * [(H - H_obs)/H0 + (2*tau/H0) * dH/dt]
This commit adds an extra term when inverting for C_c, for the following reason:
Basal shear stress for a Coulomb law is proportional to N * C_c, where N can be small when H ~ Hf.
We want to relax (N * C_c) such that H approaches a steady value without oscillating.
Thus the prognostic equation should be:
1/(N*C_c) * d(N*C_c)/dt = -(1/tau) * [(H - H_obs)/H0 + (2*tau/H0) * dH/dt]
Using the product rule on the LHS gives a term of the form (1/N)(dN/dt).
Move this term to the RHS and set dN/dt = dN/dh * dh/dt.
With N = (rhoi*g*H) * (1 - Hf/H)^p, we can show dN/dh = N * [(1 - p)/H + p/(H - Hf)].
So we add a term of this form on the RHS of the prognostic equation.
The result is to increase the rate of change of C_c near the grounding line
when ice is thinning and retreating. Ideally, this should help with stability
by reducing the chance of overshooting the GL. In practice, it may only delay
the onset of instability, if the thickness target is an unstable state (as may
be the case in the Southern Weddell Sea Embayment).
I also increased the maximum value of flow_factor_basin from 3.0 to 5.0.
In several basins, the flow factor reaches the new maximum as the ice thickens,
suggesting that we might need more melting in these basins.
This commit is answer-changing when inverting for flow_factor_basin or coulomb_c.
The new Antarctic input datasets contain fill values to indicate missing data for several fields. The fill value is a large positive number, currently 9.969213+36. To avoid errors, we need to check for fill values at initialization and replace with harmless values. This commit includes a new interface, check_fill_values, in module glide_setup. The subroutines in this interface (one for 2D fields and one for 3D fields) replace the standard netCDF fill value (declared in module glimmer_paramets) with 0.0. Optionally, the user can specify a different fill value and replacement value. This commit includes calls to check_fill_values for bheatflx, smb, and thermal_forcing, depending on the config options. With these calls, the codes runs with new input files without breaking. Note: If running with the old Glide dycore and an input file with fill values, we would need to add similar calls in glide_initialise.
There are two optional fields, f_effecpress_bwat and f_effecpress_ocean_p, that can be used to make the effective pressure N adjust gradually instead of abruptly. These fields are typically initialized to 1.0 (corresponding to no reduction in N), and then reduced during the run, decreasing N. The previous logic was to set the fields to 1 on any run that isn't a restart. However, this does not work for a run that is formally an initial run (restart = 0) but is launched from a restart file. In this case, the fields are read correctly from the restart file, but then erroneously overwritten. The fix is to check whether either field has a nonzero value in one or more cells. If so, then the field must have been read from the input file, and it is *not* overwritten. With this fix, a forward run initialized from the restart time slice at the end of a spin-up behaves correctly.
For Antarctic spin-ups, we have usually inverted for deltaT_basin, a basin-wide
temperature correction factor. Often, this leads to large thickness biases of
opposite sign within ice shelves. Part of the shelf is too thick, part is too thin,
and the temperture correction tries to compromise, but with mixed results.
This commit adds a new option, which_ho_deltaT_ocn, to apply a 2D temperature
correction factor. That is, the factor can have a different value in each floating
grid cell. Option 0, the default, does nothing; option 1 inverts for deltaT_ocn;
and option 2 applies the inversion value from an earlier spin-up.
The inversion is done in module glissade_inversion, in a new subroutine called
glissade_inversion_deltatT_ocn. This subroutine is called at the same point in the code
(during the diagnostic solve) as the bmlt_basin inversion scheme. The two options
are mutually exclusive.
The inversion works as follows. The change in the temperature correction Tc is given by
dTc/dt = (1/tau) * [T0*(H - Hobs)/H0 - (Tc - Tc_rlx)]
where T0 is a temperature scale, tau is a timescale, H0 is a thickness scale,
and Tc_rlx is a 2D field toward which Tc relaxes. For now, Tc_rlx = 0 everywhere.
Optionally, H0 is replaced by max(Hobs,H0), so the adjustment is slower where Hobs is large.
The first term raises Tc where H > Hobs and lowers Tc where H < Hobs.
The second term raises Tc where Tc < Tc_rlx and lowers Tc where Tc > Tc_rlx.
Thus, Tc can come into balance between the thickness error term and the relaxation tendency.
It does not keep going up or down indefinitely when there is a nonzero thickness error.
Note that there is no term proportional to dH/dt.
By default, tau = deltaT_ocn_timescale = 100 yr
H0 = deltaT_ocn_thck_scale = 100 m
T0 = deltaT_ocn_temp_scale = 2 deg C
These are new config parameters.
To make the Tc field less noisy, the thickness target field is smoothed, by default,
with a Laplacian smoother.
Previously, there was a 2D field called deltaT_basin. I renamed this as deltaT_ocn,
since now the ocean temperature correction does not have to be uniform within a
given basin. There is a new variable called deltaT_basin_avg, with dimension(nbasin).
This is the basin-average value used for the ISMIP6 nonlocal melt scheme.
Now it is fairly common to have bmlt_float < 0 (i.e., freeze-on).
This can happen wherever the local correction is negative enough to give
negative local thermal forcing, while the basin average thermal forcing
is still positive.
I also made some changes in basal friction inversion (C_c or C_p inversion):
* I removed the dvelo_sfc_dt term that didn't work well.
* I added logical variables fixed_thck_scale and fixed_velo_scale.
When true, we use fixed parameters in the denominators to set the size of the tendency terms.
When false, we use the local observed values to set the size of the tendency terms.
The default is true, which maintains the older convention.
I changed the calving logic to prevent the mask from moving a full cell beyond the
observed calving front at initialization. A few months ago, we added some logic
so that in regions where the observed ice speed is nonzero, the calving mask
is set to 0 (i.e., no calving). However, the mask was set to zero in ice-free ocean cells
with nonzero speed at *any* of 4 neighboring vertices. Now, the mask is set to zero
only in ice-free ocean cells with nonzero speed at all 4 neighboring vertices.
This removes a band at the CF where the ice is much thicker than the observed H = 0,
driving large, unphysical local melt rates.
…erms
This commit adds a new option, which_ho_flow_enhancement_factor, to support inversion in
each ice-covered grid cell for a flow enhancement factor, E, to match a thickness target.
E is a scalar of O(1) that multiplies the temperature-dependent flow factor A.
Most Antarctic runs to date have set E = 1 for grounded ice, with a smaller value
(e.g., 0.4 or 0.5) for floating ice. Higher E corresponds to softer ice and faster flow.
I removed an older inversion option, which_ho_flow_factor_basin, which inverted
for a large-scale flow enhancement factor (E) in each ocean basin.
There are three choices for which_ho_flow_enhancement_factor:
* 0 = uniform value of E for all grounded ice, and a uniform value for all floating ice.
These two values are often but not always different. Each has a default value of 1.0.
* 1 = invert for E(i,j) in ice-covered cells
* 2 = apply E(i,j) from an external file, usually based on a previous inversion
For option 1, each cell is initialized with one of the uniform parameters, depending on
whether it is grounded or floating. During the inversion, E is adjusted.
Where H > Hobs, E increases, and where H < Hobs, E decreases.
The inversion is carried out in subroutine glissade_inversion_flow_enhancement_factor.
The logic is similar to the subroutines that invert for Cc or Cp, and for dT_ocn.
The config parameter 'flow_factor' is now called 'flow_factor_ground'.
In any config file that sets flow_factor in the [parameters] section, the user should
change the name, although this is critical only for runs with flow_factor_ground /= 1.
The parameter 'flow_factor_float' still has the same name.
Until now, we have typically used one inversion method for grounded ice (e.g., Cc inversion),
and another for floating ice (e.g., inverting for dT_ocn either locally or at the basin scale).
With this commit, it is possible to invert simulataneously for E, Cc and dT_ocn, applying
two methods (or three, for partly grounded cells) in the same cell.
We mediate between the different inversions by adding a relaxation term to each inversion equation.
For the inversions based on local thickness (Cc or Cp, dT_ocn, and E) there are now three terms on the RHS:
(1) a term proportional to (H - Hobs)/(H0*tau), where H0 and tau are config parameters,
(2) a term proportional to (dH/dt)/H0, and
(3) a term that nudges the inverted field back toward a relaxation target, to discourage
excursions too far away from physically reasonable values.
For E, the relaxation targets are flow_factor_ground for grounded ice and
flow_factor_float for floating ice. This allows the inversion to reach a steady state
where dE/dt and dH/dt are small, giving a balance between term 1 (proportional to
the thickness error) and term 3 (the relaxation term). The greater the thickness error,
the greater the departure of E from its relaxation target.
I added a new option, which_ho_coulomb_c_relax, to set the relaxation target for Cc.
The target can be uniform for all cells (option 1), or a function of bed elevation (option 2).
For option 2, the defaults are 0.40 at high bed elevations and 0.10 at low elevations.
For inversion of dT_ocn, I added a dH/dt term to be consistent with the other options.
In new Antarctic spin-ups, the results look good. The thickness biases and GL errors
are smaller then in previous spin-ups. The relaxation terms prevent clumping
of inversion parameters at extreme values, e.g. coulomb_c = 1 or 0.001.
For these runs I used the following config parameters:
* coulomb_c inversion:
babc_thck_scale = 100 m
babc_timescale = 500 yr
babc_relax_factor = 0.05
* flow_enhancement_factor inversion:
flow_enhancement_thck_scale = 100 m
flow_enhancement_timescale = 500 yr
flow_enhancement_relax_factor = 0.5
* deltaT_ocean inversion:
deltaT_ocn_thck_scale = 100 m
deltaT_ocn_timescale = 100 yr
deltaT_ocn_temp_scale = 2.0 degC
The relaxation factors (or T0, for the case of deltaT_ocn inversion)
control the relative size of the relaxation term compared to the first two terms.
This commit changes 'flow_factor' to 'flow_factor_ground' in config files for several test cases, consistent with the name change in the previous commit.
Several months ago, I added an option to set calving_mask = 0 (i.e., do not calve) in ice-free ocean cells with a nonzero initial velocity. The nonzero velocity was taken to indicate that the cell was ice-covered at some time in recent years. This commit comments out that code. Now, the default behavior for mask-based calving is to calve ice in all cells that are initially ice-free ocean. It is still possible for the user to specify a calving mask in the input file.
This commit adds a new option, which_ho_deltaT_ocn = HO_DELTAT_OCN_DTHCK_DT = 3, to set deltaT_ocn based on dthck_dt_obs, the observed rate of ice thickness change. This option would typically be used after a long spin-up. It is assumed that by the end of the spin-up, the ice thickness is close to equilibrium. The adjustments to deltaT_ocn transform the ice from equilibrium to a state that captures the observed present-day thinning of floating ice. To make these adjustments, the user does a diagnostic restart starting from the end of the spin-up. Here, 'diagnostic' means that the final time, tend, is the same as the time of the restart file. Thus, CISM computes a new deltaT_ocn field without stepping the model forward. If not present in the restart file from the spin-up, the dthck_dt_obs field should be added. Also, in the config file, the user should set forcewrite_restart = T in the [options] section. The forcewrite_restart option is new. The default is false; setting it to true forces CISM to write a time slice to the restart file after updating deltaT_ocn. The time slice with the new deltaT_ocn field can then be used to initiate a forward run. To set deltaT_ocn, I added a subroutine, ismip6_set_deltaT_ocn, in module glissade_bmlt_float. When which_ho_deltaT_ocn = 3, this subroutine is called during the diagnostic solve. For any of the three ISMIP6 sub-shelf melt schemes, CISM increases deltaT_ocn in all cells where dthck_dt_obs < 0. This warming increases the melt rate in floating grid cells by an amount that matches dthck_dt_obs. For the nonlocal and nonlocal-slope schemes, we iterate to convergence, updating the basin-average thermal forcing on each iteration. I tested the new scheme with Tim's dH/dt (= dthck_dt_obs) dataset derived from Smith et al. (2020). I verified that the increased melt rates match the observed dthck_dt for all three ISMIP6 basal melt parameterizations (local, nonlocal, and nonlocal-slope). I reduced deltaT_ocn_maxval from 10 C to 5 C. I also changed the diagnostics for calving vs. basal melting in calving-front (CF) cells. In runs with a calving mask and deltaT_ocn inversion, CISM often diagnoses a large melt rate in CF cells. For diagnostics, it is better to think of this melt as part of the calving. With this commit, any basal melting in CF cells is used to increment the calving-rate diagnostic and then is zeroed out in the output. As a result, up to several hundred Gt/yr of basal melting can be reclassified as calving.
This commit includes two changes to Antarctic runs:
(1) Where floating ice is thinner than observed and deltaT_ocn is being driven negative,
cap deltaT_ocn to prevent the net thermal forcing from falling below zero.
This prevents the thermal forcing from remaining below zero if and when deltaT_ocn
becomes more positive. Answers will change as a result.
(2) Added a new logical config option, enable_acab_dthck_dt_correction.
The default is false. When this option is true, we add (-dthck_dt_obs) to the SMB
where the ice is floating and dthck_dt_obs < 0. This means, in many cases,
that a nonzero melt rate is needed to counteract the corrected SMB during spin-up
with inversion for deltaT_ocn.
During the subsequent forward run, we turn off inversion and set this option back to false.
Then the nonzero melt rate will drive melting and thinning similar to dthck_dt_obs.
This is similar to the method Tim has been using. It is an alternative to setting
deltaT_ocn based on dthck_dt_obs after the ice sheet has already been spun up.
…ojections I added a 2D field called dthck_dt_obs_basin. For each cell in a basin, this field contains the basin-average value of dthck_dt_obs for that basin, where the average is computed beneath ice that is floating in observations. When running with option enable_acab_dthck_dt_correction = T, the applied correction is now equal to the basin-average value of (-dthck_dt_obs), instead of the local value. The correction is zeroed out in basins with dthck_dt_obs > 0. The new field is written to the restart file when inverting for deltaT_ocn (which is always the case if enable_acab_dthck_dt_correction = T). Alos, I added a config option called deltaT_ocn_extrapolate, which is F by default. When set to T we compute the basin-average value of deltaT_ocn beneath floating ice, and then extrapolate that value to grounded ice in the basin. Should the grounding line retreat, the newly floating ice sees a basin-average value instead of zero. For the ISMIP6 Antarctic spin-ups, I set this option to T during inversion for deltaT_ocn (which_ho_deltaT_ocn = 1). The resulting values are written to the restart file. When the model is run forward using the values obtained during inversion (which_ho_deltaT_ocn = 2), the option is set back to F. (The code does this automatically if the user forgets.) The basin-average values of deltaT_ocn are read from the restart file at the start of the forward run, and these values stay the same throughout the run. The motivation is as follows. The sub-ice-shelf melt rate is sensitive to gamma0. As gamma0 increases, there is more melting. During inversion, however, the ice thickness is constrained by observations. If higher gamma0 leads to thinner ice, the inversion will compensate by computing negative values of deltaT_ocn. When the model is run forward, we have to decide how to set deltaT_ocn for grounded ice that becomes afloat. The simplest thing is to set deltaT_ocn = 0, implying that we apply the observation-based thermal forcing without a correction. In this case, higher gamma0 will lead to faster, perhaps unrealistic retreat. By setting deltaT_ocn to the basin-average values, we are correcting the thermal forcing in an appropriate way, based on the hypothesis that gamma0 is realistic but the observations may have biases. Higher gamma0 implies more negative values of deltaT_ocn, acting as a negative feedback on retreat.
This commit adds a config parameter called f_ground_threshold, part of the calving derived type. This parameter replaces a hard-wired parameter of the same name that was used in the remove_icebergs subroutine, and another called isthmus_f_ground_threshold in the remove_isthmuses subroutine. The user can now set this threshold in the config file. The default is f_ground_threshold = 0.10, consistent with the old hard-wired value for icebergs. A higher threshold means that when locating possible icebergs, a higher value of f_ground_cell (f_ground_cell >= f_ground_threshold) is needed to seed the fill. For isthmus removal, weakly grounded cells (f_ground < f_ground_threshold) and not just floating cells can be identified as isthmuses. Option 2 for force_retreat (remove floating ice marked by a mask) has changed. Instead of removing only floating ice, this options now removes ice with f_ground_cell < f_ground_threshold, so that weakly grounded cells and not just floating cells can be removed. This change improves stability by preventing small islands of weakly grounded cells from surviving when nearbly floating ice is removed. In subroutine glissade_grounded_fraction, we no longer set f_ground_cell = 1 in all land-based cells. Instead, f_ground_cell is based on the value in each quadrant, as for marine-based cells. This is answer-changing. Together, these changes improve stability when applying ice_fraction_retreat_mask, the shelf-collapse mask read in for ISMIP6 Antarctic experiments. For the ISMIP6 Antarctic 2300 projections, we used the previous commit for spin-up and for forward runs without shelf collapse. This commit was used for forward runs with shelf collapse.
This commit comes after a rebase. The code to which I rebased includes 'model_id' arguments when adding restart variables. In the course of the rebase, I missed a few. This commit puts them back so the code will build.
|
The code itself provides BFB results with main for the tested that were run. However, due to modification in the configuration options related to this PR, modifications needed to be done in the configuration files of some tests highlighted below.
Also, the README files and runCISM scripts will need to be reviewed as they were tailored to Cheyenne. But this is minor and does not impact the code or the run of the tests themselves. These changes can be addressed later. |
The basal_physics branch has some modified config options, which means that the mismip+.config.template file needs to be changed in order to run the MISMIP+ experiments correctly. The diffs are: > diff mismip+.config.template.orig mismip+.config.template < which_ho_effecpress = 3 # 3 = ocean connection < coulomb_c = 0.5 < powerlaw_c = 1.e4 # constant for the power law stress --- > coulomb_c_const = 0.5 > powerlaw_c_const = 1.e4 # constant for the power law stress That is, the uniform values of coulomb_c and powerlaw_c need 'const' in the variable names, and which_ho_effecpress changes from 3 to the default value of 0. On the basal_physics branch, setting p_ocean_penetration to a nonzero value (p = 1 for MISMIP+) changes N near the grounding line regardless of the value of which_ho_effecpress. I also made some minor changes in the MISMIP+ README file, and turned off the verbose_inversion and verbose_calving options. I verified that MISMIP+ results (including the Spinup, Ice0, and Ice1r experiments) are BFB compared to the main branch. In the log files, the calving and basal melt fluxes are now partitioned differently for Ice1r, but this is a diagnostic change that doesn't affect any state variables. I left several possible MISMIP+ changes for later: * Use the restart = 2 option for the experiments that follow the Spinup. This option is not yet in main but will come in with an upcoming merge. * Modify the setup script so as to put a working Derecho run script in each directory. Will need to load some modules (like my derecho-intel-modules script). * Add another script that launches each run script?
When I compared the basal_physics branch to the main branch using the new LIVVkit, the stream test failed. This is because which_ho_babc = 2 no longer has the same meaning. In the stream config file, I set which_ho_babc = 15 to turn on the till yield stress option used in the stream test. With this change, the stream test passes. Also turned off verbose options for calving and inversion.
|
@gunterl and @Katetc, I've finished testing this branch and comparing to main. All the standard LIVV tests pass. I confirmed that MISMIP+ results are BFB, given the config file changes described above by @gunterl. These changes are now checked in. I ran some additional DIVA dome tests and found that the agreement with main is not BFB but is within machine roundoff, which I'd say is good enough. Would one of you like to do the merge? |
|
@Katetc, Thank you! |
TimvandenAkker
pushed a commit
that referenced
this pull request
Jan 22, 2025
Basal hydrology scheme, new basal friction and inversion options
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This PR includes three main lines of effort:
(1) updates to the slab test
(2) a steady-state flux routing basal hydrology scheme
(3) various basal friction and inversion options to support Antarctic ice sheet experiments
The slab test is as described in Section 3.4 of Robinson et al. (TC, 2022). It was used in that paper
to test the stability of various solvers (SIA, SSA, DIVA, L1L2) in the simple case of an infinite slab
of ice flowing down an inclined plane. Various test parameters can be specified when launching
the run script. See the README file for details.
The hydrology scheme is based on the idea that all water input at the bed is conservatively routed down
the hydraulic gradient until it reaches the ice margin. Steady state means that for each grid cell, the water output is equal to the input. Water can be routed to one downstream neighbor (D8), two neighbors (Dinf), or all downstream neighbors (FD8). The scheme uses the efficient algorithm of Planchon & Darboux (2001) to fill depressions. It is fully parallel, unlike the similar scheme in old Glimmer. The scheme has been tested for the Greenland and Antarctic ice sheets; more testing and tuning will follow.
The basal friction and inversion changes include the following:
Results from this branch were compared to main for several standard test cases.
Using LIVVkit (with several dome, ISMIP-HOM, shelf, and stream tests), all tests pass. The stream test requires a config file change.
MISMIP+ results are also BFB, with several config file changes.
DIVA dome tests are not quite BFB, but agree within machine roundoff. I think this is acceptable, given the large number of changes on this branch.