This workflow is based on using the python package
phydrus which was developed for HYDRUS1D
v4.08 and is not compatible with the last official release of HYDRUS1D
v4.17.0140. In order to overcome this issue the R package kwb.hydrus1d is
currently developed by Kompetenzzentrum Wasser Berlin gGmbH, which will
be compatible with the last offical HYDRUS1D
v4.17.0140 release.
The workflow below is thus out-of-date and only kept for historical purposes. For the current workflow using HYDRUS1D v4.17.0140 please checkout the Modelling with R vignette.
Install R Package
# Enable this universe
options(repos = c(
kwbr = 'https://kwb-r.r-universe.dev',
CRAN = 'https://cloud.r-project.org'
))
# Install R package
install.packages('flextreat.hydrus1d')Install Phython Package phydrus
is_windows <- Sys.info()[1] == "Windows"
library(flextreat.hydrus1d)
### Compile Hydrus1D v4.08 for windows
### based on: https://github.com/phydrus/source_code/tree/main/source
if (is_windows) {
exe_path <- flextreat.hydrus1d::create_hydrus_exe()
writeLines(exe_path, "exe_path.txt")
}
env <- "flextreat"
kwb.python::conda_py_install(
env_name = env,
pkgs = list(
conda = c("python=3.9", "openpyxl"),
py = "https://github.com/phydrus/phydrus/zipball/dev/79823e797662b9838dda13888e3ead962c6699d0"
)
)
reticulate::use_miniconda(condaenv = env, required = TRUE)
tdir <- tempdir()
fs::dir_create(tdir)
env_yml <- kwb.python::conda_export(condaenv = env, export_dir = tdir)
file.copy(env_yml, basename(env_yml), overwrite = TRUE)
# For reproducing the Python environment used for the `Hydrus1D` modelling you can
# check the `env_yml` file. Keep in mind that `phydrus` was developed
# for the last open-source release (`hydrus v4.0.8`) and is not compatible with the
# latest `hydrus v4.17.0140`
# (see: https://github.com/phydrus/source_code/issues/1#issuecomment-1034675697).
env_yml
data(atmosphere)
write.csv(atmosphere, file = "atmosphere.csv", row.names = FALSE)# Import the necessary modules
import os
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import phydrus as ps
###############################################################################
# 0 General Model Settings
###############################################################################
exe_path = open("exe_path.txt").read().splitlines()
exe_path = ''.join(exe_path)
exe_path
ml = ps.model.Model(
exe_name = exe_path,
ws_name = "test",
name = "Test",
description = "Test for FlexTreat",
length_unit = "cm",
time_unit = "days",
mass_units = "mmol",
print_screen = True
)
ml.add_time_info(
tinit = 0, # Initial time of the simulation [T].
tmax = 2, # Final time of the simulation [T].
print_times = True
# Set to True. if information of pressure head, water contents,
# temperatures, and concentrations in observation nodes, and
# the water and solute fluxes is to be printed at a constant
# time interval of 1 time unit.
)
ml.add_waterflow(
top_bc = 3, #Upper Boundary Condition:
#0 = Constant Pressure Head.
#1 = Constant Flux.
#2 = Atmospheric Boundary Condition with Surface Layer.
#3 = Atmospheric Boundary Condition with Surface Run Off.
#4 = Variable Pressure Head.
#5 = Variable Pressure Head/Flux.
bot_bc = 4 #Lower Boundary Condition:
#0 = Constant Pressure Head.
#1 = Constant Flux.
#2 = Variable Pressure Head.
#3 = Variable Flux.
#4 = Free Drainage.
#5 = Deep Drainage.
#6 = Seepage Face.
#7 = Horizontal Drains.
)
## Add the process solute transport
# ml.add_solute_transport(
# lupw = True
# )
# Add materials
m = ml.get_empty_material_df(n=2)
m.loc[0:2] = [
[0.0449, 0.3833, 0.0385, 3.2783, 619.18, 0.5],
[0.04985, 0.377475, 0.03485, 3.8696, 1188.05, 0.5]
]
ml.add_material(m)
profile = ps.create_profile(
top = 0, # float, optional
# Top of the soil column.
bot = [-60, -200], # float or list of float, optional
# Bottom of the soil column. If a list is provided, multiple layers are
# created and other arguments need to be of the same length (e.g. mat).
dx = 1, # float, optional
# Size of each grid cell
h = -500, # float, optional
# Initial values of the pressure head.
lay = [1,2], # int or list of int, optional
# Subregion number (for mass balance calculations).
mat = [1,2], # int or list of int, optional
# Material number (for heterogeneity).
conc = 1e-10 # float, option
)
ml.add_profile(profile)
ml.add_obs_nodes(
depths = [-30, -60, -90, -120, -160, -200] # list of ints
# List of floats denoting the depth of the nodes. The depth is
# defined in the same length unit as selected in ml.model function.
# The function defines the closest node to the desired depth.
)
# sol = ml.get_empty_solute_df()
# sol.loc[:, "beta"] = 1.0
# ml.add_solute(
# data = sol,
# difw = 1e-9, #float, optional
# #Ionic or molecular diffusion coefficient in free water, Dw [L2T-1].
# difg = 0 #float, optional
# #Ionic or molecular diffusion coefficient in gas phase, Dg [L2T-1]
# )
atm = pd.read_csv("atmosphere.csv", index_col=0)
atm["cBot"] = 0.0
atm["cTop"] = 0.0
atm.loc[1, "cTop"] = 1.0
atm[0:2]
ml.add_atmospheric_bc(
atmosphere = atm[0:2], #pandas.DataFrame
#Pandas DataFrame with at least the following columns: tAtm,
#Prec, rSoil, rRoot, hCritA, rB, hB, hT, tTop, tBot, and Ampl.
ldailyvar = False, #bool, optional
#True if HYDRUS-1D is to generate daily variations in evaporation
#and transpiration. False otherwise.
lsinusvar = False, #bool, optional
#True if HYDRUS-1D is to generate sinusoidal variations in
#precipitation. False otherwise.
llai = False, # bool, optional
# Boolean indicating that potential evapotranspiration is
# to be divided into potential evaporation and potential
# transpiration using eq. (2.75) of the manual.
rextinct = 0.463, # float, optional
# A constant for the radiation extinction by the canopy
# (rExtinct=0.463) [-]. only used when lLai is True.
hcrits = 0, # float, optional
# Maximum allowed pressure head at the soil surface [L]. Default is
# 1e+30.
tatm = 0, # float, optional
# Time for which the i-th data record is provided [T].
hcrita = 100000.0, # float, optional
# Absolute value of the minimum allowed pressure head at the soil
# surface [L].
rb = 0, # float, optional
# Bottom flux [LT-1] (set equal to 0 if KodBot is positive, or if
# one of the logical variables qGWLF, FreeD or SeepF is .true.).
hb = 0, # float, optional
# Groundwater level [L], or any other prescribed pressure head
# boundary condition as indicated by a positive value of KodBot
# (set equal to 0 if KodBot is negative, or if one of the logical
# variables qGWLF, FreeD or SeepF is .true.).
ht = 0, # float, optional
# Prescribed pressure head [L] at the surface (set equal to 0 if
# KodBot is negative).
ttop = 0, # float, optional
# Soil surface temperature [oC] (omit if both lTemp and lChem are
# equal to .false.).
tbot = 0, # float, optional
# Soil temperature at the bottom of the soil profile [oC] (omit if
# both lTemp and lChem are equal to .false., set equal to zero if
# kBotT=0).
ampl = 0 # float, optional
# Temperature amplitude at the soil surface [K] (omit if both lTemp
# and lChem are equal to .false.).
)
# Write the input and check if the model simulates
ml.write_input()