Added a steady-state basal water routing scheme #38
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
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 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.
whlipscomb
force-pushed
the
lipscomb/basal_hydrology
branch
from
0378fe2
to
d56ece0
Compare
|
@Katetc, this PR can be removed. The basal hydrology changes, like the slab test changes, were brought to main as part of basal_physics4. |
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.
These commits add a steady-state flux routing hydrology scheme. It routes basal water conservatively down the hydraulic gradient from source cells (where bmlt > 0) to the ice sheet boundary. It is fully parallel, adapted from a similar serial scheme that was written some years ago for Glimmer. By conserving basal water and allowing the water to flow through cells with basal temperature below the pressure melting point, this scheme improves on the current local till model used in many Greenland simulations. To use the new scheme with Glissade, the user sets 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 main computational difficulties arise from interior depressions that can be sinks for meltwater. These regions are handled using the algorithm of Planchon and Darboux (2001). In depressions, the value of the hydraulic potential is raised just enough to give every cell a downhill route to the ice sheet boundary, and no farther. This algorithm parallelizes fairly efficiently.
The scheme computes the steady-state water flux through each cell. Given this flux and an expression for water velocity, it is easy to compute basal water depth. However, these depths are small (~ 1 mm or less). It is unclear whether effective pressure should depend on the flux or the nominal depth.
From a given cell, water can be routed downhill to one down-gradient neighbor, two neighbors, or up to eight neighbors depending on the value of ho_flux_routing_scheme, a new config option. More testing is needed to determine which routing scheme is most realistic.
The scheme has been tested by Sarah Bradley and colleagues in LGM runs with a North American Ice Sheet, and shown to give rise to ice streams where none were present before. It will next be tested for Greenland and Antarctica.
Answers are BFB, apart from runs that use the new hydrology scheme