From a6665a89a7a1ffbfa4e19a63678632d9fbee9037 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sun, 1 Dec 2024 11:10:00 +1100 Subject: [PATCH] Include pure Julia implementation of IHACRES and rename node type BilinearNode -> IHACRESBilinearNode Remove level calculation which can be done elsewhere. --- src/Nodes/IHACRES/IHACRESExpuhNode.jl | 44 +-- src/Nodes/IHACRES/IHACRESNode.jl | 446 +++++++++++------------ src/Nodes/IHACRES/components/cmd.jl | 153 ++++++++ src/Nodes/IHACRES/components/flow.jl | 202 ++++++++++ src/Nodes/IHACRES/components/p_and_et.jl | 112 ++++++ 5 files changed, 712 insertions(+), 245 deletions(-) create mode 100644 src/Nodes/IHACRES/components/cmd.jl create mode 100644 src/Nodes/IHACRES/components/flow.jl create mode 100644 src/Nodes/IHACRES/components/p_and_et.jl diff --git a/src/Nodes/IHACRES/IHACRESExpuhNode.jl b/src/Nodes/IHACRES/IHACRESExpuhNode.jl index 50578e3..cbfbea6 100644 --- a/src/Nodes/IHACRES/IHACRESExpuhNode.jl +++ b/src/Nodes/IHACRES/IHACRESExpuhNode.jl @@ -17,17 +17,17 @@ Base.@kwdef mutable struct ExpuhNode{P, A<:AbstractFloat} <: IHACRESNode storage_coef::P = Param(2.9, bounds=(0.2, 10.0)) - level_params::Array{P, 1} = [ - Param(-0.01, bounds=(-10.0, -0.01)), # p1 - Param(0.8, bounds=(0.0, 1.5)), # p2 - Param(4.5, bounds=(0.0, 20.0)), # p3 - Param(5.0, bounds=(1.0, 10.0)), # p4 - Param(0.35, bounds=(0.0, 1.0)), # p5 - Param(1.41, bounds=(-2.0, 2.0)), # p6 - Param(-1.45, bounds=(-2.5, 0.0)), # p7 - Param(6.75, bounds=(0.0, 10.0)), # p8 - Param(150.0, bounds=(50.0, 200.0)) # ctf - ] + # level_params::Array{P, 1} = [ + # Param(-0.01, bounds=(-10.0, -0.01)), # p1 + # Param(0.8, bounds=(0.0, 1.5)), # p2 + # Param(4.5, bounds=(0.0, 20.0)), # p3 + # Param(5.0, bounds=(1.0, 10.0)), # p4 + # Param(0.35, bounds=(0.0, 1.0)), # p5 + # Param(1.41, bounds=(-2.0, 2.0)), # p6 + # Param(-1.45, bounds=(-2.5, 0.0)), # p7 + # Param(6.75, bounds=(0.0, 10.0)), # p8 + # Param(150.0, bounds=(50.0, 200.0)) # ctf + # ] storage::Array{A} = [100.0] quick_store::Array{A} = [0.0] @@ -47,17 +47,17 @@ function ExpuhNode(name::String, spec::Dict) node_params = spec["parameters"] n_lparams = n.level_params s_lparams = spec["level_params"] - node_params["level_params"] = Param[ - Param(s_lparams[1], bounds=n_lparams[1].bounds) - Param(s_lparams[2], bounds=n_lparams[2].bounds) - Param(s_lparams[3], bounds=n_lparams[3].bounds) - Param(s_lparams[4], bounds=n_lparams[4].bounds) - Param(s_lparams[5], bounds=n_lparams[5].bounds) - Param(s_lparams[6], bounds=n_lparams[6].bounds) - Param(s_lparams[7], bounds=n_lparams[7].bounds) - Param(s_lparams[8], bounds=n_lparams[8].bounds) - Param(s_lparams[9], bounds=n_lparams[9].bounds) - ] + # node_params["level_params"] = Param[ + # Param(s_lparams[1], bounds=n_lparams[1].bounds) + # Param(s_lparams[2], bounds=n_lparams[2].bounds) + # Param(s_lparams[3], bounds=n_lparams[3].bounds) + # Param(s_lparams[4], bounds=n_lparams[4].bounds) + # Param(s_lparams[5], bounds=n_lparams[5].bounds) + # Param(s_lparams[6], bounds=n_lparams[6].bounds) + # Param(s_lparams[7], bounds=n_lparams[7].bounds) + # Param(s_lparams[8], bounds=n_lparams[8].bounds) + # Param(s_lparams[9], bounds=n_lparams[9].bounds) + # ] node_params["storage"] = [node_params["initial_storage"]] delete!(node_params, "initial_storage") diff --git a/src/Nodes/IHACRES/IHACRESNode.jl b/src/Nodes/IHACRES/IHACRESNode.jl index d1cee59..5af421e 100644 --- a/src/Nodes/IHACRES/IHACRESNode.jl +++ b/src/Nodes/IHACRES/IHACRESNode.jl @@ -1,11 +1,14 @@ using Parameters using ModelParameters +include("./components/p_and_et.jl") +include("./components/cmd.jl") +include("./components/flow.jl") abstract type IHACRESNode <: NetworkNode end -Base.@kwdef mutable struct BilinearNode{P, A<:AbstractFloat} <: IHACRESNode +Base.@kwdef mutable struct IHACRESBilinearNode{P, A<:AbstractFloat} <: IHACRESNode const name::String const area::A @@ -24,17 +27,17 @@ Base.@kwdef mutable struct BilinearNode{P, A<:AbstractFloat} <: IHACRESNode storage_coef::P = Param(2.9, bounds=(1e-10, 10.0)) alpha::P = Param(0.1, bounds=(1e-5, 1 - 1/10^9)) - const level_params::Array{P, 1} = [ - Param(-0.01, bounds=(-10.0, -0.01)), # p1 - Param(0.8, bounds=(0.0, 1.5)), # p2 - Param(4.5, bounds=(0.0, 20.0)), # p3 - Param(5.0, bounds=(1.0, 10.0)), # p4 - Param(0.35, bounds=(0.0, 1.0)), # p5 - Param(1.41, bounds=(-2.0, 2.0)), # p6 - Param(-1.45, bounds=(-2.5, 0.0)), # p7 - Param(6.75, bounds=(0.0, 10.0)), # p8 - Param(150.0, bounds=(50.0, 200.0)) # ctf - ] + # const level_params::Array{P, 1} = [ + # Param(-0.01, bounds=(-10.0, -0.01)), # p1 + # Param(0.8, bounds=(0.0, 1.5)), # p2 + # Param(4.5, bounds=(0.0, 20.0)), # p3 + # Param(5.0, bounds=(1.0, 10.0)), # p4 + # Param(0.35, bounds=(0.0, 1.0)), # p5 + # Param(1.41, bounds=(-2.0, 2.0)), # p6 + # Param(-1.45, bounds=(-2.5, 0.0)), # p7 + # Param(6.75, bounds=(0.0, 10.0)), # p8 + # Param(150.0, bounds=(50.0, 200.0)) # CTF (height in local datum where stream will cease to flow) + # ] storage::Array{A} = [100.0] # CMD quick_store::Array{A} = [0.0] @@ -43,11 +46,8 @@ Base.@kwdef mutable struct BilinearNode{P, A<:AbstractFloat} <: IHACRESNode effective_rainfall::Array{A} = [] et::Array{A} = [] inflow::Array{A} = [] - level::Array{A} = [] + # level::Array{A} = [] gw_store::Array{A} = [0.0] - - # cache arrays - cache::NamedTuple{(:i_cache, :r_cache), <:Tuple{Vector{Float64}, Vector{Float64}}} = (i_cache=zeros(3), r_cache=zeros(2)) end @@ -65,7 +65,7 @@ function prep_state!(node::IHACRESNode, timesteps::Int64)::Nothing node.effective_rainfall = fill(0.0, timesteps) node.et = fill(0.0, timesteps) node.inflow = fill(0.0, timesteps) - node.level = fill(0.0, timesteps) + # node.level = fill(0.0, timesteps) resize!(node.gw_store, timesteps+1) node.gw_store[2:end] .= 0.0 @@ -74,25 +74,25 @@ function prep_state!(node::IHACRESNode, timesteps::Int64)::Nothing end -function BilinearNode(name::String, spec::Dict) - n = create_node(BilinearNode, name, spec["area"]) +function IHACRESBilinearNode(name::String, spec::Dict) + n = create_node(IHACRESBilinearNode, name, spec["area"]) node_params = copy(spec["parameters"]) - if haskey(spec, "level_params") - n_lparams = n.level_params - s_lparams = spec["level_params"] - node_params["level_params"] = Param[ - Param(s_lparams[1], bounds=n_lparams[1].bounds) - Param(s_lparams[2], bounds=n_lparams[2].bounds) - Param(s_lparams[3], bounds=n_lparams[3].bounds) - Param(s_lparams[4], bounds=n_lparams[4].bounds) - Param(s_lparams[5], bounds=n_lparams[5].bounds) - Param(s_lparams[6], bounds=n_lparams[6].bounds) - Param(s_lparams[7], bounds=n_lparams[7].bounds) - Param(s_lparams[8], bounds=n_lparams[8].bounds) - Param(s_lparams[9], bounds=n_lparams[9].bounds) - ] - end + # if haskey(spec, "level_params") + # n_lparams = n.level_params + # s_lparams = spec["level_params"] + # node_params["level_params"] = Param[ + # Param(s_lparams[1], bounds=n_lparams[1].bounds) + # Param(s_lparams[2], bounds=n_lparams[2].bounds) + # Param(s_lparams[3], bounds=n_lparams[3].bounds) + # Param(s_lparams[4], bounds=n_lparams[4].bounds) + # Param(s_lparams[5], bounds=n_lparams[5].bounds) + # Param(s_lparams[6], bounds=n_lparams[6].bounds) + # Param(s_lparams[7], bounds=n_lparams[7].bounds) + # Param(s_lparams[8], bounds=n_lparams[8].bounds) + # Param(s_lparams[9], bounds=n_lparams[9].bounds) + # ] + # end node_params["storage"] = [node_params["initial_storage"]] delete!(node_params, "initial_storage") @@ -121,16 +121,16 @@ end """ - BilinearNode(name::String, area::Float64, d::Float64, d2::Float64, e::Float64, f::Float64, + IHACRESBilinearNode(name::String, area::Float64, d::Float64, d2::Float64, e::Float64, f::Float64, a::Float64, b::Float64, s_coef::Float64, alpha::Float64, store::Float64, quick::Float64, slow::Float64, gw_store::Float64) Create a IHACRES node that adopts the bilinear CMD module. """ -function BilinearNode(name::String, area::Float64, d::Float64, d2::Float64, e::Float64, f::Float64, +function IHACRESBilinearNode(name::String, area::Float64, d::Float64, d2::Float64, e::Float64, f::Float64, a::Float64, b::Float64, s_coef::Float64, alpha::Float64, store::Float64, quick::Float64, slow::Float64, gw_store::Float64) - n = create_node(BilinearNode, name, area) + n = create_node(IHACRESBilinearNode, name, area) update_params!(n, d, d2, e, f, a, b, s_coef, alpha) @@ -143,9 +143,16 @@ function BilinearNode(name::String, area::Float64, d::Float64, d2::Float64, e::F end -function update_state!(s_node::IHACRESNode, storage::Float64, e_rainfall::Float64, et::Float64, - qflow_store::Float64, sflow_store::Float64, outflow::Float64, - level::Float64, gw_store::Float64)::Nothing +function update_state!( + s_node::IHACRESNode, + storage::F, + e_rainfall::F, + et::F, + qflow_store::F, + sflow_store::F, + outflow::F, + gw_store::F +)::Nothing where {F<:Float64} push!(s_node.storage, storage) push!(s_node.effective_rainfall, e_rainfall) push!(s_node.et, et) @@ -153,7 +160,7 @@ function update_state!(s_node::IHACRESNode, storage::Float64, e_rainfall::Float6 push!(s_node.quick_store, qflow_store) push!(s_node.slow_store, sflow_store) - push!(s_node.level, level) + # push!(s_node.level, level) push!(s_node.gw_store, gw_store) return nothing @@ -168,7 +175,6 @@ function update_state!( sflow_store::Float64, inflow::Float64, outflow::Float64, - level::Float64, gw_store::Float64 )::Nothing s_node.storage[ts+1] = storage @@ -181,14 +187,16 @@ function update_state!( s_node.effective_rainfall[ts] = e_rainfall s_node.et[ts] = et - s_node.level[ts] = level + # s_node.level[ts] = level return nothing end """ - run_node!(node::IHACRESNode, climate::Climate, timestep::Int; - inflow=nothing, extraction=nothing, exchange=nothing) + run_node!( + node::IHACRESNode, climate::Climate, ts::Int; + inflow=nothing, extraction=nothing, exchange=nothing + ) Run a specific node for a specified time step. @@ -200,8 +208,10 @@ Run a specific node for a specified time step. - `extraction::DataFrame` : Time series of water orders (expects column of `_releases`) - `exchange::DataFrame` : Time series of groundwater flux """ -function run_node!(node::IHACRESNode, climate::Climate, ts::Int; - inflow=nothing, extraction=nothing, exchange=nothing) +function run_node!( + node::IHACRESNode, climate::Climate, ts::Int; + inflow=nothing, extraction=nothing, exchange=nothing +) if checkbounds(Bool, node.outflow, ts) if node.outflow[ts] != undef # already ran for this time step so no need to run @@ -220,11 +230,11 @@ end """ run_timestep!( - node::BilinearNode, climate::Climate, timestep::Int, + node::IHACRESBilinearNode, climate::Climate, timestep::Int, inflow::Float64, extraction::Float64, exchange::Float64 ) -Run the given BilinearNode for a time step. +Run the given IHACRESBilinearNode for a time step. """ function run_timestep!( node::IHACRESNode, climate::Climate, timestep::Int; @@ -235,7 +245,7 @@ function run_timestep!( return run_timestep!(node, ts, P, ET, in_flow, ext, flux) end -function run_timestep!(node, P, ET, ts; inflow=nothing, extraction=nothing, exchange=nothing) +function run_timestep!(node::IHACRESNode, P, ET, ts; inflow=nothing, extraction=nothing, exchange=nothing) in_flow = timestep_value(ts, node.name, "inflow", inflow) ext = timestep_value(ts, node.name, "extraction", extraction) flux = timestep_value(ts, node.name, "exchange", exchange) @@ -251,94 +261,84 @@ function run_ihacres!( inflow::F, ext::F, gw_exchange::F -)::Tuple{F, F} where {F<:Float64} +)::F where {F<:Float64} current_store::F = s_node.storage[ts] quick_store::F = s_node.quick_store[ts] slow_store::F = s_node.slow_store[ts] gw_store::F = s_node.gw_store[ts] - cache_res = s_node.cache.i_cache - @ccall IHACRES.calc_ft_interim_cmd( - cache_res::Ptr{Cdouble}, - current_store::Cdouble, - rain::Cdouble, - s_node.d.val::Cdouble, - s_node.d2.val::Cdouble, - s_node.alpha.val::Cdouble - )::Cvoid - (mf, e_rainfall, recharge) = cache_res - - et::F = @ccall IHACRES.calc_ET( - s_node.e.val::Cdouble, - evap::Cdouble, - mf::Cdouble, - s_node.f.val::Cdouble, - s_node.d.val::Cdouble - )::Cdouble - - cmd::F = @ccall IHACRES.calc_cmd( - current_store::Cdouble, - rain::Cdouble, - et::Cdouble, - e_rainfall::Cdouble, - recharge::Cdouble - )::Cdouble - - @ccall IHACRES.calc_ft_flows( - cache_res::Ptr{Cdouble}, - quick_store::Cdouble, - slow_store::Cdouble, - e_rainfall::Cdouble, - recharge::Cdouble, - s_node.area::Cdouble, - s_node.a.val::Cdouble, - s_node.b.val::Cdouble - )::Cvoid - (nq_store, ns_store, outflow) = cache_res - - routing_res = s_node.cache.r_cache - @ccall IHACRES.routing( - routing_res::Ptr{Cdouble}, - gw_store::Cdouble, - s_node.storage_coef.val::Cdouble, - inflow::Cdouble, - outflow::Cdouble, - ext::Cdouble, - gw_exchange::Cdouble - )::Cvoid - - (gw_store, outflow) = routing_res - - level_params = Vector{Float64}(s_node.level_params) - level::Float64 = @ccall IHACRES.calc_ft_level( - outflow::Cdouble, - level_params::Ptr{Cdouble} - )::Cdouble - - update_state!(s_node, ts, cmd, e_rainfall, et, nq_store, ns_store, inflow, outflow, level, gw_store) + (mf, e_rainfall, recharge) = calc_ft_interim_cmd( + current_store, + rain, + s_node.d.val::F, + s_node.d2.val::F, + s_node.alpha.val::F + ) - return outflow, level + et::F = calc_ET_from_E( + s_node.e.val::F, + evap, + mf, + s_node.f.val::F, + s_node.d.val::F + ) + + cmd::F = calc_cmd( + current_store, + rain, + et, + e_rainfall, + recharge + ) + + (nq_store, ns_store, outflow) = calc_ft_flows( + quick_store, + slow_store, + e_rainfall, + recharge, + s_node.area, + s_node.a.val::F, + s_node.b.val::F + ) + + (gw_store, outflow) = routing( + gw_store, + s_node.storage_coef.val::F, + inflow, + outflow, + ext, + gw_exchange + ) + + # level::Float64 = calc_ft_level( + # outflow, + # Vector{Float64}(s_node.level_params) + # ) + + update_state!(s_node, ts, cmd, e_rainfall, et, nq_store, ns_store, inflow, outflow, gw_store) + + return outflow end """ - param_info(node::IHACRESNode; with_level::Bool = true)::Tuple + param_info(node::IHACRESNode)::Tuple Extract node parameter names, values, and bounds for IHACRESNode types. """ -function param_info(node::IHACRESNode; with_level::Bool = false)::Tuple +function param_info(node::IHACRESNode)::Tuple tmp = Model(node) values = collect(tmp[:val]) bounds = collect(tmp[:bounds]) param_names = collect(tmp[:fieldname]) - if with_level - level_param_vals = map(x->x.val, node.level_params) - append!(values, level_param_vals) + # if with_level + # level_param_vals = map(x->x.val, node.level_params) + # append!(values, level_param_vals) - level_param_bounds = map(x->x.bounds, node.level_params) - append!(bounds, level_param_bounds) - end + # level_param_bounds = map(x->x.bounds, node.level_params) + # append!(bounds, level_param_bounds) + # end return param_names, values, bounds end @@ -357,7 +357,7 @@ function run_node_with_temp!(sn::StreamfallNetwork, nid::Int64, climate::Climate end -function run_node_with_temp!(node::BilinearNode, climate::Climate; +function run_node_with_temp!(node::IHACRESBilinearNode, climate::Climate; inflow=nothing, extraction=nothing, exchange=nothing) timesteps = sim_length(climate) @inbounds for ts in 1:timesteps @@ -369,7 +369,7 @@ function run_node_with_temp!(node::BilinearNode, climate::Climate; end -function run_node_with_temp!(node::BilinearNode, climate::Climate, timestep::Int64; +function run_node_with_temp!(node::IHACRESBilinearNode, climate::Climate, timestep::Int64; inflow=nothing, extraction=nothing, exchange=nothing) ts = timestep P, T = climate_values(node, climate, ts) @@ -381,7 +381,7 @@ end """ - run_node_with_temp!(s_node::BilinearNode, + run_node_with_temp!(s_node::IHACRESBilinearNode, rain::Float64, temp::Float64, inflow::Float64, @@ -394,7 +394,7 @@ end Run node with temperature data to calculate outflow and update state. """ -function run_node_with_temp!(s_node::BilinearNode, +function run_node_with_temp!(s_node::IHACRESBilinearNode, rain::Float64, temp::Float64, inflow::Float64, @@ -409,82 +409,82 @@ function run_node_with_temp!(s_node::BilinearNode, slow_store = s_node.slow_store[end] gw_store = s_node.gw_store[end] - interim_results = [0.0, 0.0, 0.0] - @ccall IHACRES.calc_ft_interim_cmd(interim_results::Ptr{Cdouble}, - current_store::Cdouble, - rain::Cdouble, - s_node.d::Cdouble, - s_node.d2::Cdouble, - s_node.alpha::Cdouble)::Cvoid - - @assert any(isnan.(interim_results)) == false - (mf, e_rainfall, recharge) = interim_results - - et::Float64 = @ccall IHACRES.calc_ET_from_T( - s_node.e::Cdouble, - temp::Cdouble, - mf::Cdouble, - s_node.f::Cdouble, - s_node.d::Cdouble - )::Cdouble - - cmd::Float64 = @ccall IHACRES.calc_cmd( - mf::Cdouble, - rain::Cdouble, - et::Cdouble, - e_rainfall::Cdouble, - recharge::Cdouble - )::Cdouble - - flow_results = [0.0, 0.0, 0.0] - @ccall IHACRES.calc_ft_flows( - flow_results::Ptr{Cdouble}, - quick_store::Cdouble, - slow_store::Cdouble, - e_rainfall::Cdouble, - recharge::Cdouble, - s_node.area::Cdouble, - s_node.a::Cdouble, - s_node.b::Cdouble - )::Cvoid - - @assert any(isnan.(flow_results)) == false - (nq_store, ns_store, outflow) = flow_results - - routing_res = [0.0, 0.0] - @ccall IHACRES.routing( - routing_res::Ptr{Cdouble}, - gw_store::Cdouble, - s_node.storage_coef::Cdouble, - inflow::Cdouble, - outflow::Cdouble, - ext::Cdouble, - gw_exchange::Cdouble)::Cvoid - - @assert any(isnan.(routing_res)) == false - (gw_store, outflow) = routing_res - - level_params = Array{Float64}(s_node.level_params) - level::Float64 = @ccall IHACRES.calc_ft_level( - outflow::Cdouble, - level_params::Ptr{Cdouble} - )::Cdouble + (mf, e_rainfall, recharge) = calc_ft_interim_cmd( + interim_results, + current_store, + rain, + s_node.d, + s_node.d2, + s_node.alpha + ) + + et::Float64 = calc_ET_from_T( + s_node.e, + temp, + mf, + s_node.f, + s_node.d + ) + + cmd::Float64 = calc_cmd( + mf, + rain, + et, + e_rainfall, + recharge + ) + + (nq_store, ns_store, outflow) = calc_ft_flows( + quick_store, + slow_store, + e_rainfall, + recharge, + s_node.area, + s_node.a, + s_node.b + ) + + @assert any(isnan.((nq_store, ns_store, outflow))) == false + + (gw_store, outflow) = routing( + gw_store, + s_node.storage_coef, + inflow, + outflow, + ext, + gw_exchange + ) + + # level_params = Array{Float64}(s_node.level_params) + # level::Float64 = calc_ft_level( + # outflow, + # level_params + # ) push!(s_node.inflow, inflow) - update_state(s_node, cmd, e_rainfall, et, nq_store, ns_store, outflow, level, gw_store) + update_state(s_node, cmd, e_rainfall, et, nq_store, ns_store, outflow, gw_store) return outflow, level end """ - update_params!(node::BilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64, + update_params!(node::IHACRESBilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64, a::Float64, b::Float64, s_coef::Float64, alpha::Float64)::Nothing Update model parameters. """ -function update_params!(node::BilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64, - a::Float64, b::Float64, s_coef::Float64, alpha::Float64)::Nothing +function update_params!( + node::IHACRESBilinearNode, + d::F, + d2::F, + e::F, + f::F, + a::F, + b::F, + s_coef::F, + alpha::F +)::Nothing where {F<:Float64} node.d = Param(d, bounds=node.d.bounds::Tuple) node.d2 = Param(d2, bounds=node.d2.bounds::Tuple) node.e = Param(e, bounds=node.e.bounds::Tuple) @@ -498,41 +498,41 @@ function update_params!(node::BilinearNode, d::Float64, d2::Float64, e::Float64, end -""" - update_params!(node::BilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64, - a::Float64, b::Float64, s_coef::Float64, alpha::Float64, - p1::Float64, p2::Float64, p3::Float64, p4::Float64, p5::Float64, p6::Float64, p7::Float64, p8::Float64, CTF::Float64)::Nothing - -Update all parameters. -""" -function update_params!(node::BilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64, - a::Float64, b::Float64, s_coef::Float64, alpha::Float64, - p1::Float64, p2::Float64, p3::Float64, p4::Float64, p5::Float64, - p6::Float64, p7::Float64, p8::Float64, CTF::Float64)::Nothing - node.d = Param(d, bounds=node.d.bounds::Tuple) - node.d2 = Param(d2, bounds=node.d2.bounds::Tuple) - node.e = Param(e, bounds=node.e.bounds::Tuple) - node.f = Param(f, bounds=node.f.bounds::Tuple) - node.a = Param(a, bounds=node.a.bounds::Tuple) - node.b = Param(b, bounds=node.b.bounds::Tuple) - node.storage_coef = Param(s_coef, bounds=node.storage_coef.bounds::Tuple) - node.alpha = Param(alpha, bounds=node.alpha.bounds::Tuple) - - n_lparams = node.level_params - node.level_params = [ - Param(p1, bounds=n_lparams[1].bounds::Tuple) - Param(p2, bounds=n_lparams[2].bounds::Tuple) - Param(p3, bounds=n_lparams[3].bounds::Tuple) - Param(p4, bounds=n_lparams[4].bounds::Tuple) - Param(p5, bounds=n_lparams[5].bounds::Tuple) - Param(p6, bounds=n_lparams[6].bounds::Tuple) - Param(p7, bounds=n_lparams[7].bounds::Tuple) - Param(p8, bounds=n_lparams[8].bounds::Tuple) - Param(CTF, bounds=n_lparams[9].bounds::Tuple) - ] - - return nothing -end +# """ +# update_params!(node::IHACRESBilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64, +# a::Float64, b::Float64, s_coef::Float64, alpha::Float64)::Nothing + +# Update all parameters. +# """ +# function update_params!(node::BilinearNode, d::F, d2::F, e::F, f::F, +# a::F, b::F, s_coef::F, alpha::F, +# p1::F, p2::F, p3::F, p4::F, p5::F, +# p6::F, p7::F, p8::F, CTF::F +# )::Nothing where {F<:Float64} +# node.d = Param(d, bounds=node.d.bounds::Tuple) +# node.d2 = Param(d2, bounds=node.d2.bounds::Tuple) +# node.e = Param(e, bounds=node.e.bounds::Tuple) +# node.f = Param(f, bounds=node.f.bounds::Tuple) +# node.a = Param(a, bounds=node.a.bounds::Tuple) +# node.b = Param(b, bounds=node.b.bounds::Tuple) +# node.storage_coef = Param(s_coef, bounds=node.storage_coef.bounds::Tuple) +# node.alpha = Param(alpha, bounds=node.alpha.bounds::Tuple) + +# n_lparams = node.level_params +# node.level_params = [ +# Param(p1, bounds=n_lparams[1].bounds::Tuple) +# Param(p2, bounds=n_lparams[2].bounds::Tuple) +# Param(p3, bounds=n_lparams[3].bounds::Tuple) +# Param(p4, bounds=n_lparams[4].bounds::Tuple) +# Param(p5, bounds=n_lparams[5].bounds::Tuple) +# Param(p6, bounds=n_lparams[6].bounds::Tuple) +# Param(p7, bounds=n_lparams[7].bounds::Tuple) +# Param(p8, bounds=n_lparams[8].bounds::Tuple) +# Param(CTF, bounds=n_lparams[9].bounds::Tuple) +# ] + +# return nothing +# end # TODO # update_bounds!() @@ -549,7 +549,7 @@ function reset!(s_node::IHACRESNode)::Nothing s_node.slow_store = [s_node.slow_store[1]] s_node.gw_store = [s_node.gw_store[1]] s_node.outflow = [] - s_node.level = [] + # s_node.level = [] s_node.effective_rainfall = [] s_node.et = [] s_node.inflow = [] diff --git a/src/Nodes/IHACRES/components/cmd.jl b/src/Nodes/IHACRES/components/cmd.jl new file mode 100644 index 0000000..c21597a --- /dev/null +++ b/src/Nodes/IHACRES/components/cmd.jl @@ -0,0 +1,153 @@ +""" + calc_cmd(prev_cmd::Float64, rainfall::Float64, et::Float64, effective_rainfall::Float64, recharge::Float64)::Float64 + +Calculate Catchment Moisture Deficit (CMD). + +Min value of CMD is 0.0 and is represented in mm depth. +A value of 0 indicates that the catchment is fully saturated. +A value greater than 0 means that there is a moisture deficit. +""" +function calc_cmd( + prev_cmd::Float64, rainfall::Float64, et::Float64, + effective_rainfall::Float64, recharge::Float64 +)::Float64 + cmd = prev_cmd + et + effective_rainfall + recharge - rainfall # units in mm + + return max(0.0, cmd) +end + +""" + calc_linear_interim_cmd(cmd::Float64, param_d::Float64, rainfall::Float64)::Float64 + +Calculate interim CMD (M_f) in its linear form. + +Based on HydroMad implementation and details in references. + +# Arguments +- `cmd`: previous Catchment Moisture Deficit (M_{k-1}) +- `param_d`: model parameter factor `d` +- `rainfall`: rainfall for current time step in mm + +# Returns +interim CMD (M_f) + +# References +- Croke, B.F.W., Jakeman, A.J. 2004. A catchment moisture deficit module for the IHACRES rainfall-runoff model, + Environmental Modelling & Software, 19(1), pp. 1–5. doi: 10.1016/j.envsoft.2003.09.001 +- Croke, B.F.W., Jakeman, A.J. 2005. Corrigendum to "A Catchment Moisture Deficit module for the IHACRES + rainfall-runoff model [Environ. Model. Softw. 19 (1) (2004) 1–5]" Environmental Modelling & Software, + 20(7), p. 997. doi: 10.1016/j.envsoft.2004.11.004 +""" +function calc_linear_interim_cmd(cmd::Float64, param_d::Float64, rainfall::Float64)::Float64 + if cmd < param_d + return cmd * exp(-rainfall / param_d) + elseif cmd < (param_d + rainfall) + return param_d * exp((-rainfall + cmd - param_d) / param_d) + end + + return cmd - rainfall +end + +""" + calc_trig_interim_cmd(cmd::Float64, param_d::Float64, rainfall::Float64)::Float64 + +Calculate interim CMD (M_f) in its trigonometric form. + +Based on HydroMad implementation and details in references. + +# Arguments +- `cmd`: previous Catchment Moisture Deficit (M_{k-1}) +- `param_d`: model parameter `d` +- `rainfall`: rainfall for current time step in mm + +# Returns +- interim CMD (M_f) +""" +function calc_trig_interim_cmd(cmd::Float64, param_d::Float64, rainfall::Float64)::Float64 + if cmd < param_d + Mf = 1.0 / tan((cmd / param_d) * (π / 2.0)) + return (2.0 * param_d / π) * atan(1.0 / (π * rainfall / (2.0 * param_d) + Mf)) + elseif rainfall < (param_d + rainfall) + return (2.0 * param_d / π) * atan(2.0 * param_d / (π * (param_d - cmd + rainfall))) + else + return cmd - rainfall + end +end + +""" + calc_ft_interim_cmd(cmd::Float64, rain::Float64, d::Float64, d2::Float64, alpha::Float64) + +Direct port of original Fortran implementation to calculate interim CMD (M_f). +Calculates estimates of effective rainfall and recharge as a by-product. + +# Arguments +- `cmd`: previous Catchment Moisture Deficit (M_{k-1}) +- `rain`: rainfall for time step in mm +- `d`: flow threshold value +- `d2`: scaling factor applied to `d` +- `alpha`: scaling factor + +# Returns +A tuple of (interim CMD value, effective rainfall, recharge) +""" +function calc_ft_interim_cmd(cmd::Float64, rain::Float64, d::Float64, d2_factor::Float64, alpha::Float64)::Tuple{Float64, Float64, Float64} + # Initialize variables + d2 = d * d2_factor + tmp_cmd = cmd + e_rain = 0.0 + recharge = 0.0 + + # Early return if no rain + if isapprox(rain, 0.0) + return (cmd, e_rain, recharge) + end + + # Calculate Mf (new CMD value) + if tmp_cmd > (d2 + rain) + # CMD never reaches d2, so all rain is effective + Mf = tmp_cmd - rain + else + # Handle more complex cases + if tmp_cmd > d2 + tmp_rain = rain - (tmp_cmd - d2) # leftover rain after reaching d2 threshold + tmp_cmd = d2 + else + tmp_rain = rain + end + + # Calculate threshold value d1a + d1a = d * (2.0 - exp(-(rain / 50.0)^2)) + + if tmp_cmd > d1a + eps = d2 / (1.0 - alpha) + + # Amount of rain necessary to get to threshold d + depth_to_d = eps * log((alpha + tmp_cmd / eps) / (alpha + d1a / eps)) + + if depth_to_d >= tmp_rain + lam = exp(tmp_rain * (1.0 - alpha) / d2) # lambda + epsilon = alpha * eps + + Mf = tmp_cmd / lam - epsilon * (1.0 - 1.0 / lam) + e_rain = 0.0 + else + if tmp_cmd > d1a + tmp_rain = tmp_rain - depth_to_d + end + + tmp_cmd = d1a + gamma = (alpha * d2 + (1.0 - alpha) * d1a) / (d1a * d2) + Mf = tmp_cmd * exp(-tmp_rain * gamma) + e_rain = alpha * (tmp_rain + 1.0 / d1a / gamma * (Mf - tmp_cmd)) + end + else + gamma = (alpha * d2 + (1.0 - alpha) * d1a) / (d1a * d2) + Mf = tmp_cmd * exp(-tmp_rain * gamma) + e_rain = alpha * (tmp_rain + 1.0 / d1a / gamma * (Mf - tmp_cmd)) + end + + recharge = rain - (cmd - Mf) - e_rain + end + + return (Mf, e_rain, recharge) +end diff --git a/src/Nodes/IHACRES/components/flow.jl b/src/Nodes/IHACRES/components/flow.jl new file mode 100644 index 0000000..d736721 --- /dev/null +++ b/src/Nodes/IHACRES/components/flow.jl @@ -0,0 +1,202 @@ +approx_zero(x) = x + one(x) ≈ one(x) + +""" + convert_taus(tau::Float64, v::Float64) + +Convert τ to α and β + +# Arguments +- `tau` : time constant τ +- `v` : proportion of flow v (between 0 and 1) + +# Returns +Tuple of (alpha, beta) +""" +@inline function convert_taus(tau::Float64, v::Float64)::Tuple{Float64,Float64} + return (exp(-1.0 / tau), v * (1.0 - alpha)) +end + +""" + calc_stores(tau::Float64, prev_store::Float64, v::Float64, e_rainfall::Float64) + +Calculate flow stores. +""" +function calc_stores(tau::Float64, prev_store::Float64, v::Float64, e_rainfall::Float64)::Float64 + alpha, beta = convert_taus(tau, v) + store = (beta * e_rainfall) + alpha * prev_store + return max(0.0, store) +end + +""" + calc_flows(prev_quick::Float64, prev_slow::Float64, v_s::Float64, + e_rainfall::Float64, area::Float64, tau_q::Float64, tau_s::Float64) + +Calculate quick and slow flow, and outflow using a linear routing module. + +Assumes components are in parallel. + +# Arguments +- `prev_quick` : previous quick flow in ML/day +- `prev_slow` : previous slow flow in ML/day +- `v_s` : proportional amount that goes to slow flow. v_s <= 1.0 +- `e_rainfall` : effective rainfall for t +- `area` : catchment area +- `tau_q` : Time constant, quick flow τ value +- `tau_s` : Time constant slow flow τ value + +# Returns +Tuple of (quickflow, slowflow, outflow) in ML/day + +# References +- https://wiki.ewater.org.au/display/SD45/IHACRES-CMD+-+SRG +- Croke, B.F.W., Jakeman, A.J. 2004 + A catchment moisture deficit module for the IHACRES rainfall-runoff model, + Environmental Modelling & Software, 19(1), pp. 1–5. + doi: 10.1016/j.envsoft.2003.09.001 +- Croke, B.F.W., Jakeman, A.J. 2005 + Corrigendum to "A Catchment Moisture Deficit module for the IHACRES + rainfall-runoff model [Environ. Model. Softw. 19 (1) (2004) 1–5]" + Environmental Modelling & Software, 20(7), p. 997. + doi: 10.1016/j.envsoft.2004.11.004 +""" +function calc_flows(prev_quick::Float64, prev_slow::Float64, v_s::Float64, + e_rainfall::Float64, area::Float64, tau_q::Float64, tau_s::Float64 + )::Tuple{Float64,Float64,Float64} + v_q = 1.0 - v_s # proportional quick flow + areal_rainfall = e_rainfall * area + quick = calc_stores(tau_q, prev_quick, v_q, areal_rainfall) + slow = calc_stores(tau_s, prev_slow, v_s, areal_rainfall) + outflow = quick + slow + return (quick, slow, outflow) +end + +""" + routing(gw_vol::Float64, storage_coef::Float64, inflow::Float64, flow::Float64, + irrig_ext::Float64, gw_exchange::Float64=0.0) + +Stream routing, taking into account groundwater interactions and water extractions. + +# Arguments +- gw_vol : groundwater store at t-1 +- storage_coef : groundwater storage factor +- inflow : incoming streamflow (flow from previous node) +- flow : outflow for the node (local flow) +- irrig_ext : volume of irrigation extraction in ML +- gw_exchange : groundwater flux. Defaults to 0.0 + Negative values represent infiltration into aquifer + +# Returns +Tuple of (gw_store, streamflow) in ML/day +""" +function routing(gw_vol::Float64, storage_coef::Float64, inflow::Float64, flow::Float64, + irrig_ext::Float64, gw_exchange::Float64=0.0)::Tuple{Float64,Float64} + tmp_gw_store = gw_vol + (inflow + flow + gw_exchange) - irrig_ext + + if tmp_gw_store > 0.0 + # Account for interaction with groundwater system + c1 = exp(-storage_coef) + gw_store = c1 * tmp_gw_store + outflow = (1.0 - c1) * tmp_gw_store + else + # Groundwater level is below stream, so no baseflow occurs + gw_store = tmp_gw_store + outflow = 0.0 + end + + return (gw_store, outflow) +end + +""" + calc_outflow(flow::Float64, extractions::Float64) + +Simple routing modifier to account for extractions. + +# Arguments +- `flow` : sum of quickflow and slowflow in ML/day +- `extractions` : water extractions that occurred in ML/day + +# Returns +(flow - extractions), minimum possible is 0.0 +""" +function calc_outflow(flow::Float64, extractions::Float64)::Float64 + return max(0.0, flow - extractions) +end + +""" + calc_ft_flows(prev_quick::Float64, prev_slow::Float64, e_rain::Float64, + recharge::Float64, area::Float64, a::Float64, b::Float64) + +Unit Hydrograph module ported from Fortran. + +# Arguments +- `prev_quick` : previous quickflow storage +- `prev_slow` : previous slowflow storage +- `e_rain` : effective rainfall in mm +- `recharge` : recharge amount in mm +- `area` : catchment area in km² +- `a` : quickflow storage coefficient, inverse of τ_q such that a := (1/τ_q) +- `b` : slowflow storage coefficient, inverse of τ_s such that b := (1/τ_s) + +# Returns +Tuple of (quick_store, slow_store, outflow) +""" +function calc_ft_flows(prev_quick::Float64, prev_slow::Float64, e_rain::Float64, + recharge::Float64, area::Float64, a::Float64, b::Float64 + )::Tuple{Float64,Float64,Float64} + tmp_calc = max(0.0, prev_quick + (e_rain * area)) + + if tmp_calc > 0.0 + alpha = exp(-a) + beta = (1.0 - alpha) * tmp_calc + quick_store = alpha * tmp_calc + outflow = beta + else + quick_store = tmp_calc + outflow = 0.0 + end + + @assert outflow >= 0.0 "Calculating quick store: Outflow cannot be negative: $outflow" + + slow_store = prev_slow + (recharge * area) + if slow_store > 0.0 + alpha = exp(-b) + beta = (1.0 - alpha) * slow_store + outflow += beta + slow_store = alpha * slow_store + end + + @assert outflow >= 0.0 "Calculating slow store: Outflow cannot be negative: $outflow" + + return (quick_store, slow_store, outflow) +end + +""" + calc_ft_level(outflow::Float64, level_params::Vector{Float64}) + +Calculate stream level from outflow using level parameters. + +# Arguments +- `outflow` : stream outflow +- `level_params` : array of 9 parameters [p1, p2, p3, p4, p5, p6, p7, p8, CTF] + +# Returns +Calculated stream level +""" +function calc_ft_level(outflow::Float64, level_params::Vector{Float64})::Float64 + p1, p2, p3, p4, p5, p6, p7, p8, CTF = level_params + + if approx_zero(outflow) + return CTF # Return cease to flow level + end + + level = (exp(p1) * outflow^p2) * + (1.0 / ( + 1.0 + ((outflow / p3)^p4)^(p5/p4) * exp(p6 / (1.0 + exp(-p7*p8)) * outflow^p7) + )) + level = max(level, 0.0) + level += CTF # add Cease to Flow (base height of stream in local datum) + + @assert level >= 0.0 "Stream level cannot be below 0.0, got $level" + + return level +end diff --git a/src/Nodes/IHACRES/components/p_and_et.jl b/src/Nodes/IHACRES/components/p_and_et.jl new file mode 100644 index 0000000..25526a2 --- /dev/null +++ b/src/Nodes/IHACRES/components/p_and_et.jl @@ -0,0 +1,112 @@ +export calc_effective_rainfall, calc_ET_from_E, calc_ET, calc_ET_from_T + +""" + calc_effective_rainfall(rainfall::Float64, cmd::Float64, d::Float64, d2::Float64, n::Float64=0.1) + +Estimate effective rainfall. + +# References +- Croke, B.F.W., Jakeman, A.J. 2004 + A catchment moisture deficit module for the IHACRES rainfall-runoff model, + Environmental Modelling & Software, 19(1), pp. 1–5. + doi: 10.1016/j.envsoft.2003.09.001 + +- Croke, B.F.W., Jakeman, A.J. 2005 + Corrigendum to "A Catchment Moisture Deficit module for the IHACRES + rainfall-runoff model [Environ. Model. Softw. 19 (1) (2004) 1–5]" + Environmental Modelling & Software, 20(7), p. 997. + doi: https://doi.org/10.1016/j.envsoft.2004.11.004 + +# Arguments +- `rainfall`: rainfall for time step +- `cmd`: previous CMD value +- `d`: threshold value +- `d2`: scaling factor applied to `d` +- `n`: scaling factor (default = 0.1). Default value is suitable for most cases (Croke & Jakeman, 2004) + +# Returns +- Effective rainfall value +""" +function calc_effective_rainfall(rainfall::Float64, cmd::Float64, d::Float64, + d2::Float64, n::Float64=0.1)::Float64 + scaled_d = d * d2 + + if cmd > scaled_d + e_rainfall = rainfall + else + f1 = min(1.0, cmd / d) + f2 = min(1.0, cmd / scaled_d) + e_rainfall = rainfall * ((1.0 - n) * (1.0 - f1) + (n * (1.0 - f2))) + end + + return max(0.0, e_rainfall) +end + +""" + calc_ET_from_E(e::Float64, evap::Float64, Mf::Float64, f::Float64, d::Float64) + +Calculate evapotranspiration from evaporation. + +# Arguments +- `e`: temperature to PET conversion factor (a stress threshold) +- `evap`: evaporation for given time step +- `Mf`: Catchment Moisture Deficit prior to accounting for ET losses (`Mf`) +- `f`: calibrated parameter that acts as a multiplication factor on `d` +- `d`: flow threshold factor + +# Returns +- Estimate of ET +""" +function calc_ET_from_E(e::Float64, evap::Float64, Mf::Float64, f::Float64, + d::Float64)::Float64 + param_g = f * d + et = e * evap + + if Mf > param_g + et *= min(1.0, exp((1.0 - Mf/param_g)*2.0)) + end + + return max(0.0, et) +end + +# """ +# calc_ET(e::Float64, evap::Float64, Mf::Float64, f::Float64, d::Float64) + +# !!! warning +# Deprecated function - use [`calc_ET_from_E`](@ref) instead. +# """ +# function calc_ET(e::Float64, evap::Float64, Mf::Float64, f::Float64, d::Float64)::Float64 +# return calc_ET_from_E(e, evap, Mf, f, d) +# end + +""" + calc_ET_from_T(e::Float64, T::Float64, Mf::Float64, f::Float64, d::Float64) + +Calculate evapotranspiration based on temperature data. + +Parameters `f` and `d` are used to calculate `g`, the value of the CMD +which the ET rate will begin to decline due to insufficient +water availability for plant transpiration. + +# Arguments +- `e`: temperature to PET conversion factor (a stress threshold) +- `T`: temperature in degrees C +- `Mf`: Catchment Moisture Deficit prior to accounting for ET losses (`Mf`) +- `f`: multiplication factor on `d` +- `d`: flow threshold factor + +# Returns +- Estimate of ET from temperature (for catchment area) +""" +function calc_ET_from_T(e::Float64, T::Float64, Mf::Float64, f::Float64, + d::Float64)::Float64 + # temperature can be negative, so we have a min cap of 0.0 + if T <= 0.0 + return 0.0 + end + + param_g = f * d + et = e * T * min(1.0, exp(2.0 * (1.0 - (Mf / param_g)))) + + return max(0.0, et) +end