Subset CONUS and do a ParFlow Spinup

To launch this notebook interactively in a Jupyter notebook-like browser interface, please click the “Launch Binder” button below. Note that Binder may take several minutes to launch.

Binder

This notebook will subset a HUC8 from the CONUS1 domain as in this example. However, it does not subset climate forcing data as it is an example of model initialization or spin up. So, it only runs with ParFlow. The template runscript used is conus1_pf_spinup_solid.yaml which has the necessary timing keys to run with a growth timestep and a longterm recharge mask (PmE) as opposed to a constant time step and time varying climate forcing data. The example notebook is not set to run for the entire period it would normally take to spin-up a model, you are encouraged to do your own evaluation to determine if your model is at steady state.

This notebook has two principal sections:

  1. Subset all static inputs from a CONUS1 run stored in Hydrodata

  2. Load and alter a reference run to spin-up your subset ParFlow domain

Import the required libraries

import matplotlib.pyplot as plt
import os
from parflow import Run
from parflow.tools.fs import mkdir
from parflow.tools.settings import set_working_directory
import subsettools as st
import hf_hydrodata as hf

# You need to register on https://hydrogen.princeton.edu/pin before you can use the hydrodata utilities
hf.register_api_pin("your_email", "your_pin")

1. Define variables to access datasets in Hydrodata to subset and define write paths

We will be testing with the Upper Verde watershed for this example

  • HUC: 15060202

  • Size: 6496 km^2 (ni = 112, nj = 90)

Set your variables to specify which static data you would like to subset in Hydrodata

runname = "spinup_test"

# provide a way to create a subset from the conus domain (huc list, lat/lon bbox currently supported)
hucs = ["15060202"]

# provide information about the datasets you want to access for run inputs using the data catalog
start = "2005-10-01"
wy = 2006
grid = "conus1"
run_ds = "conus1_baseline_mod"
var_ds = "conus1_domain"
P = 1
Q = 1

# set the directory paths where you want to write your subset files
home = os.path.expanduser("~")
base_dir = os.path.join(home, "subsettools_tutorial")
input_dir = os.path.join(base_dir, "inputs", f"{runname}_{grid}_{wy}WY_spinup")
output_dir = os.path.join(base_dir, "outputs")
static_write_dir = os.path.join(input_dir, "static")
mkdir(static_write_dir)
pf_out_dir = os.path.join(output_dir, f"{runname}_{grid}_{wy}WY_spinup")
mkdir(pf_out_dir)

# Set the PARFLOW_DIR path to your local installation of ParFlow.
# This is only necessary if this environment variable is not already set.
# os.environ["PARFLOW_DIR"] = "/path/to/your/parflow/installation"

# load your preferred template runscript
reference_run = st.get_template_runscript(grid, "spinup", "solid", pf_out_dir)

2. Get the desired ParFlow i/j bbox from user provided geospatial information

ij_bounds, mask = st.define_huc_domain(hucs=hucs, grid=grid)
print("ij_bounds returns [imin, jmin, imax, jmax]")
print(f"bounding box: {ij_bounds}")

nj = ij_bounds[3] - ij_bounds[1]
ni = ij_bounds[2] - ij_bounds[0]
print(f"nj: {nj}")
print(f"ni: {ni}")

ij_bounds returns [imin, jmin, imax, jmax]
bounding box: (375, 239, 487, 329)
nj: 90
ni: 112

3. Make the mask and solid file

You only do this if you providing a huc or list of hucs. Otherwise, the reference run provided is for a box domain.

mask_solid_paths = st.write_mask_solid(mask=mask, grid=grid, write_dir=static_write_dir)

Wrote mask.pfb
Wrote solidfile and mask_vtk with total z of 500 meters

4. Subset the static ParFlow inputs

Two options to subset static inputs.

  1. subset_static(): This function when provided with a variable dataset hosted on hydrodata will subset all static inputs required to do a baseline run. Note that the function will raise an error if any of the requested variables do not exist in the dataset, so we need to modify the default variable list and remove “mannings” and “flowbarrier”. Pressure is the steady state pressure.

  2. subset_press_init(): This function will write the subset pressure of the last hour in the last day before your start date in the provided timezone. If no such pressure file exists in the hydrodata run dataset specifed, no file will be written. The function assumes UTC of zero as the default, but allows the user to override.

static_paths = st.subset_static(ij_bounds, dataset=var_ds, write_dir=static_write_dir,
                                var_list=("slope_x", "slope_y", "pf_indicator", "pme", "ss_pressure_head",)
                               )

Wrote slope_x.pfb in specified directory.
Wrote slope_y.pfb in specified directory.
Wrote pf_indicator.pfb in specified directory.
Wrote pme.pfb in specified directory.
Wrote ss_pressure_head.pfb in specified directory.

press_init_filepath = st.subset_press_init(
    ij_bounds, dataset=run_ds, date=start, write_dir=static_write_dir, time_zone="UTC"
)

UTC Date: 2005-10-01 00:00:00
Wrote /home/ga6/subsettools_tutorial/inputs/spinup_test_conus1_2006WY_spinup/static/conus1_baseline_mod_2005.10.01:00.00.00_UTC0_press.pfb in specified directory.

5. Set up a baseline run from a reference yaml

This function will return the correct template yaml file to do your run based on the grid, if you’re doing spin-up and if you’re using a solid file with the necessary keys changed to spin up your subset

runscript_path = st.edit_runscript_for_subset(
    ij_bounds, runscript_path=reference_run, runname=runname
)

New runname: spinup_test provided, a new yaml file will be created
No forcing directory provided, run.Solver.CLM.MetFilePath key not set
ComputationalGrid.NY set to 90 and NX to 112
GeomInput.domaininput.InputType detected as SolidFile, no additional keys to change for subset
Updated runscript written to /home/ga6/subsettools_tutorial/outputs/spinup_test_conus1_2006WY_spinup

6. Copy over your static files to your run directory

You may only need to do this once, or you may want to copy subset static files to different run directories

st.copy_files(read_dir=static_write_dir, write_dir=pf_out_dir)

7. Change the file names in your runscript if desired

If you have changed the name of a static input file either from those used in the reference yamls provided, or have changed the name of an individual file for an ensemble or other experiment, you can change it with this function by providing the target runscript (yaml or pfidb) and the new file name(s) as an arguments. Only those arguments with a specified file name will be updated

init_press = os.path.basename(press_init_filepath)
runscript_path = st.change_filename_values(
    runscript_path=runscript_path,
    init_press=init_press,
)

Initial pressure filename changed to conus1_baseline_mod_2005.10.01:00.00.00_UTC0_press.pfb
Updated runscript written to /home/ga6/subsettools_tutorial/outputs/spinup_test_conus1_2006WY_spinup

8. Change processor topolgoy if desired and then distribute your inputs and forcings to match the new topology

runscript_path = st.dist_run(
                topo_p=P,
                topo_q=Q,
                runscript_path=runscript_path,
                dist_clim_forcing=False,
           )

No forcing directory provided, only distributing static inputs
Distributed mask.pfb with NZ 1
Distributed slope_x.pfb with NZ 1
Distributed slope_y.pfb with NZ 1
Distributed pf_indicator.pfb with NZ 5
Distributed pme.pfb with NZ 5
Distributed ss_pressure_head.pfb with NZ 5
Distributed conus1_baseline_mod_2005.10.01:00.00.00_UTC0_press.pfb with NZ 5

9. Do a baseline run.

Load in the yaml run file you’ve created which is in the same folder as your static forcings. This assumes you do not want to make any changes from the parent model you used (Ex. conus1 baseline) and will run your subset at baseline conditions. Outputs should be almost identical to the parent model at your subset location for the same time period if you make no additional changes.

set_working_directory(pf_out_dir)

# load the specified run script
run = Run.from_definition(runscript_path)
print(f"Loaded run with runname: {run.get_name()}")

# The following line shows timing info that is likely important for interacting with your spin-up run

# run.TimingInfo.BaseUnit = 1.0
# run.TimingInfo.StartCount = 0
# run.TimingInfo.StartTime = 0
# normally set to something like 1000000.0 but shorter for demo purposes
run.TimingInfo.StopTime = 48.0
run.TimingInfo.DumpInterval = 24.0

# run.TimeStep.Type = 'Growth'
# run.TimeStep.InitialStep = .01
# run.TimeStep.GrowthFactor = 1.1
# run.TimeStep.MaxStep = 250.0
# run.TimeStep.MinStep = 0.0001

run.run(working_directory=pf_out_dir)

Loaded run with runname: spinup_test

# ==============================================================================
# ParFlow directory
#  - /home/SHARED/software/parflow/3.10.0
# ParFlow version
#  - 3.10.0
# Working directory
#  - /home/ga6/subsettools_tutorial/outputs/spinup_test_conus1_2006WY_spinup
# ParFlow database
#  - spinup_test.pfidb
# ==============================================================================


# ==============================================================================
# ParFlow ran successfully 💦 💦 💦 
# ==============================================================================