Constraining subsidence to part of the domain¶
One case that has been explored in the literature with the DeltaRCM model is the case of subsidence limited to one region of the model domain [1]. This model configuration can be readily achieved with model subclassing.
Setting up the custom subclass¶
class ConstrainedSubsidenceModel(pyDeltaRCM.DeltaModel):
"""A simple subclass of DeltaModel with subsidence region constrained.
This subclass *overwrites* the `init_subsidence` method to
constrain subsidence to only one region of the model domain.
"""
def __init__(self, input_file=None, **kwargs):
# inherit base DeltaModel methods
super().__init__(input_file, **kwargs)
def init_subsidence(self):
"""Initialize subsidence pattern constrained to a tighter region.
Uses theta1 and theta2 to set the angular bounds for the
subsiding region. theta1 and theta2 are set in relation to the
inlet orientation. The inlet channel is at an angle of 0, if
theta1 is -pi/3 radians, this means that the angle to the left of
the inlet that will be included in the subsiding region is 30
degrees. theta2 defines the right angular bounds for the subsiding
region in a similar fashion.
"""
_msg = 'Initializing custom subsidence field'
self.log_info(_msg, verbosity=1)
if self._toggle_subsidence:
theta1 = -(np.pi / 3)
theta2 = 0
R1 = 0.3 * self.L # radial limits (fractions of L)
R2 = 0.8 * self.L
Rloc = np.sqrt((self.y - self.L0)**2 + (self.x - self.W / 2.)**2)
thetaloc = np.zeros((self.L, self.W))
thetaloc[self.y > self.L0 - 1] = np.arctan(
(self.x[self.y > self.L0 - 1] - self.W / 2.)
/ (self.y[self.y > self.L0 - 1] - self.L0 + 1))
self.subsidence_mask = ((R1 <= Rloc) & (Rloc <= R2) &
(theta1 <= thetaloc) &
(thetaloc <= theta2))
self.subsidence_mask[:self.L0, :] = False
self.sigma = self.subsidence_mask * self.subsidence_rate * self.dt
Now, initialize the model and look at the field. Note that the colorscale depicts the magnitude of subsidence in the model per timestep (sigma, which has units meters).
mdl = ConstrainedSubsidenceModel(toggle_subsidence=True)
fig, ax = plt.subplots()
mdl.show_attribute('sigma', grid=False)
plt.show()
Using the custom subclass with the preprocessor¶
We can configure a Preprocessor to handle a set of custom runs in conjunction with out custom pyDeltaRCM model subclass.
For example, in [1], the authors explore the impact of subsidence at various rates: 3 mm/yr, 6 mm/yr, 10 mm/yr, 25 mm/yr, 50 mm/yr, and 100 mm/yr.
We can scale these rates, assuming a model intermittency factor of 0.019, representing 7 of 365 days of flooding per year, by using the convenience function scale_relative_sea_level_rise_rate:
from pyDeltaRCM.preprocessor import scale_relative_sea_level_rise_rate
subsidence_mmyr = np.array([3, 6, 10, 25, 50, 100])
subsidence_scaled = scale_relative_sea_level_rise_rate(subsidence_mmyr, If=0.019)
Now, we use matrix expansion to set up the runs with a preprocessor. For example, in a Python script, following the definition of the subclass above, define a dictionary with a matrix key and supply to the Preprocessor:
# add a matrix with subsidence to the dict
param_dict = {}
param_dict['matrix'] = {'subsidence_rate': subsidence_scaled}
# add other configurations
param_dict.update(
{'out_dir': 'liang_2016_reproduce',
'toggle_subsidence': True,
'parallel': 3}) # we can take advantage of parallel jobs
# create the preprocessor
pp = pyDeltaRCM.Preprocessor(
param_dict,
timesteps=10000)
And finally run the jobs by specifying the model subclass as the class to use when instantiating the jobs with the preprocessor.
# run the jobs
pp.run_jobs(DeltaModel=ConstrainedSubsidenceModel)
We can check whether the runs were set up, as expected:
from matplotlib.colors import Normalize
fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(10, 4))
norm = Normalize(vmin=3, vmax=100)
for i, job in enumerate(pp.job_list):
# first convert the field to a rate
subsidence_rate_field = (job.deltamodel.sigma / job.deltamodel.dt)
# now convert to mm/yr
subsidence_rate_field = (subsidence_rate_field * 1000 *
pyDeltaRCM.shared_tools._scale_factor(If=0.019, units='years'))
# and display
im = ax.flat[i].imshow(subsidence_rate_field, norm=norm)
fig.colorbar(im, ax=ax.ravel().tolist())
plt.show()
