diff --git a/Project.toml b/Project.toml index 70f03500..044484d5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Streamfall" uuid = "09d4d480-e905-4803-959d-af438f1c2811" authors = ["ConnectedSystems"] -version = "0.6.0" +version = "0.6.1" [deps] BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209" @@ -17,6 +17,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" @@ -39,14 +40,14 @@ ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [weakdeps] GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" +GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" -GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231" [extensions] GraphMakieExt = ["Makie", "GraphMakie"] -MakieExt = "Makie" GraphPlotExt = ["StatsPlots", "GraphPlot"] +MakieExt = "Makie" PlotsExt = "StatsPlots" [compat] @@ -61,6 +62,8 @@ DataFrames = "1.6" DataStructures = "0.18.10" Dates = "1.11.0" Distributions = "0.25.100" +Glob = "1.3.1" +GraphPlot = "0.6" Graphs = "1.9" JSON = "0.20.1, 0.21" LaTeXStrings = "1.3" @@ -75,8 +78,7 @@ Setfield = "1" Statistics = "1.10, 1.11.1" StatsBase = "0.34" StatsFuns = "1.3" +StatsPlots = "0.15" YAML = "0.4" ZipFile = "0.10" julia = "1" -StatsPlots = "0.15" -GraphPlot = "0.6" diff --git a/src/Nodes/IHACRES/IHACRESBilinearNodeGW.jl b/src/Nodes/IHACRES/IHACRESBilinearNodeGW.jl new file mode 100644 index 00000000..ffbb387d --- /dev/null +++ b/src/Nodes/IHACRES/IHACRESBilinearNodeGW.jl @@ -0,0 +1,473 @@ +using Parameters +using ModelParameters + +include("./components/gw_flow.jl") + + +""" + IHACRESBilinearNodeGW + +IHACRES node with explicit groundwater storage module (IHACRES_GW). + +Extends standard IHACRES with spatially-explicit groundwater component that: +- Tracks groundwater storage volume (ML) +- Calculates baseflow from storage-discharge relationship +- Allows groundwater extraction and external fluxes +- Converts storage to water table depth for calibration to bore observations + +Based on Ivkovic et al. (2005) IHACRES_GW formulation. + +# Fields +## Identifiers +- `name::String` : Node identifier +- `area::Float64` : Catchment area (km²) - used for runoff calculations +- `aquifer_area::Float64` : Aquifer area (km²) - used for GW depth conversions (defaults to catchment area) + +## CMD Parameters (same as standard IHACRES) +- `d::P` : CMD threshold (mm), controls runoff generation +- `d2::P` : Scaling factor for second threshold +- `e::P` : ET conversion factor +- `f::P` : Plant stress threshold factor +- `alpha::P` : Effective rainfall scaling factor + +## Quick Flow Parameter +- `a::P` : Quickflow coefficient (1/τ_q) + +## Groundwater Parameters (NEW) +- `tau_s::P` : Slow flow time constant (days), controls baseflow recession +- `vs::P` : Recharge fraction (0-1), proportion of effective rainfall to GW +- `L::P` : Constant loss term (ML/day), positive=loss, negative=gain +- `specific_yield::P` : Aquifer specific yield (dimensionless, typically 0.1-0.3) + +## Elevation Parameters (for depth conversion) +- `gauge_elevation::Float64` : Stream bed elevation at gauge (mAHD) +- `bore_ground_elevation::Float64` : Ground surface elevation at bore (mAHD) + +## State Variables +- `storage::Array{Float64}` : CMD (mm), length = timesteps + 1 +- `quick_store::Array{Float64}` : Quickflow storage (ML), length = timesteps + 1 +- `gw_storage::Array{Float64}` : Groundwater storage (ML), length = timesteps + 1 +- `effective_rainfall::Array{Float64}` : Effective rainfall (mm/day), length = timesteps +- `et::Array{Float64}` : Evapotranspiration (mm/day), length = timesteps +- `quickflow::Array{Float64}` : Quickflow discharge (ML/day), length = timesteps +- `baseflow::Array{Float64}` : Baseflow discharge (ML/day), length = timesteps +- `outflow::Array{Float64}` : Total outflow (ML/day), length = timesteps +- `inflow::Array{Float64}` : Inflow from upstream (ML/day), length = timesteps + +## Calibration +- `obj_func::Function` : Objective function for calibration +""" +Base.@kwdef mutable struct IHACRESBilinearNodeGW{P,A<:AbstractFloat} <: IHACRESNode + const name::String + const area::A + aquifer_area::A = area # Defaults to catchment area if not specified + + # CMD parameters (same as standard IHACRES) + d::P = Param(200.0, bounds=(10.0, 200.0), desc="CMD threshold") + d2::P = Param(2.0, bounds=(0.0001, 10.0), desc="Scaling factor") + e::P = Param(1.0, bounds=(0.9, 1.0), desc="ET conversion factor") + f::P = Param(0.8, bounds=(0.01, 2.5), desc="Plant stress threshold") + alpha::P = Param(0.1, bounds=(0.1, 1 - 1 / 10^9), desc="Effective rainfall scaling") + + # Quick flow parameter + a::P = Param(0.9, bounds=(0.1, 10.0), desc="Quickflow coefficient (1/tau_q)") + + # Groundwater parameters (NEW) + tau_s::P = Param(10.0, bounds=(10.0, 900.0), desc="Slow flow time constant (days)") + vs::P = Param(0.5, bounds=(0.1, 0.9), desc="Recharge fraction (to GW)") + L::P = Param(0.0, bounds=(-10.0, 10.0), desc="Constant loss (ML/day)") + specific_yield::P = Param(0.15, bounds=(0.05, 0.25), desc="Aquifer specific yield") + + # Elevation parameters (fixed, not calibrated) + gauge_elevation::Float64 = 100.0 # Stream bed elevation (mAHD) + bore_ground_elevation::Float64 = 110.0 # Bore ground surface (mAHD) + + # State variables + storage::Array{A} = [100.0] # CMD (mm) + quick_store::Array{A} = [0.0] # Quickflow storage (ML) + gw_storage::Array{A} = [0.0] # Groundwater storage (ML) + effective_rainfall::Array{A} = [] # mm/day + et::Array{A} = [] # mm/day + quickflow::Array{A} = [] # ML/day + baseflow::Array{A} = [] # ML/day + outflow::Array{A} = [] # ML/day + inflow::Array{A} = [] # ML/day + + obj_func::Function = obj_func +end + + +""" + IHACRESBilinearNodeGW(name::String, spec::AbstractDict) + +Construct IHACRESBilinearNodeGW from specification dictionary (YAML/Dict). + +# Required spec fields +- `area` : Catchment area (km²) +- `aquifer_area` : Aquifer area (km²) - for depth conversions +- `initial_storage` : Initial CMD (mm) +- `parameters` : Dict with parameter values +- `gauge_elevation` : Stream bed elevation (mAHD) +- `bore_ground_elevation` : Bore ground elevation (mAHD) + +# Optional spec fields +- `initial_gw_storage` : Initial GW storage (ML), defaults to 0.0 +""" +function IHACRESBilinearNodeGW(name::String, spec::AbstractDict) + n = create_node(IHACRESBilinearNodeGW, name, spec["area"]) + node_params = copy(spec["parameters"]) + + # Set initial CMD storage + n.storage = [spec["initial_storage"]] + + # Set initial GW storage (if provided) + if haskey(spec, "initial_gw_storage") + n.gw_storage = [spec["initial_gw_storage"]] + end + + # Set aquifer area (required) + n.aquifer_area = spec["aquifer_area"] + + # Set elevation parameters + n.gauge_elevation = spec["gauge_elevation"] + n.bore_ground_elevation = spec["bore_ground_elevation"] + + # Set parameters from spec + for (k, p) in node_params + s = Symbol(k) + if p isa String + p = eval(Meta.parse(p)) + end + + try + f = getfield(n, s) + setfield!(n, s, Param(p, bounds=f.bounds)) + catch err + msg = sprint(showerror, err, catch_backtrace()) + if occursin("no field bounds", string(msg)) + setfield!(n, s, p) + else + throw(err) + end + end + end + + return n +end + + +""" + IHACRESBilinearNodeGW(name, area, d, d2, e, f, a, tau_s, vs, L, specific_yield, + gauge_elev, bore_elev, cmd_store, quick_store, gw_store) + +Construct IHACRESBilinearNodeGW with direct parameter specification. +""" +function IHACRESBilinearNodeGW( + name::String, area::Float64, + d::Float64, d2::Float64, e::Float64, f::Float64, + a::Float64, alpha::Float64, + tau_s::Float64, vs::Float64, L::Float64, specific_yield::Float64, + gauge_elev::Float64, bore_elev::Float64, + cmd_store::Float64, quick_store::Float64, gw_store::Float64 +) + n = create_node(IHACRESBilinearNodeGW, name, area) + + # Set parameters + n.d = Param(d, bounds=n.d.bounds) + n.d2 = Param(d2, bounds=n.d2.bounds) + n.e = Param(e, bounds=n.e.bounds) + n.f = Param(f, bounds=n.f.bounds) + n.a = Param(a, bounds=n.a.bounds) + n.alpha = Param(alpha, bounds=n.alpha.bounds) + n.tau_s = Param(tau_s, bounds=n.tau_s.bounds) + n.vs = Param(vs, bounds=n.vs.bounds) + n.L = Param(L, bounds=n.L.bounds) + n.specific_yield = Param(specific_yield, bounds=n.specific_yield.bounds) + + # Set elevations + n.gauge_elevation = gauge_elev + n.bore_ground_elevation = bore_elev + + # Set initial states + n.storage = [cmd_store] + n.quick_store = [quick_store] + n.gw_storage = [gw_store] + + return n +end + + +""" + prep_state!(node::IHACRESBilinearNodeGW, timesteps::Int64)::Nothing + +Prepare node for simulation by pre-allocating state arrays. +""" +function prep_state!(node::IHACRESBilinearNodeGW, timesteps::Int64)::Nothing + resize!(node.storage, timesteps + 1) + node.storage[2:end] .= 0.0 + + resize!(node.quick_store, timesteps + 1) + node.quick_store[2:end] .= 0.0 + + resize!(node.gw_storage, timesteps + 1) + node.gw_storage[2:end] .= 0.0 + + node.effective_rainfall = fill(0.0, timesteps) + node.et = fill(0.0, timesteps) + node.quickflow = fill(0.0, timesteps) + node.baseflow = fill(0.0, timesteps) + node.outflow = fill(0.0, timesteps) + node.inflow = fill(0.0, timesteps) + + return nothing +end + + +""" + reset!(node::IHACRESBilinearNodeGW)::Nothing + +Reset node state to initial conditions. +""" +function reset!(node::IHACRESBilinearNodeGW)::Nothing + node.storage = [node.storage[1]] + node.quick_store = [node.quick_store[1]] + node.gw_storage = [node.gw_storage[1]] + node.effective_rainfall = [] + node.et = [] + node.quickflow = [] + node.baseflow = [] + node.outflow = [] + node.inflow = [] + + return nothing +end + + +""" + run_timestep!(node::IHACRESBilinearNodeGW, rain::F, evap::F, ts::Int64; + inflow=nothing, extraction=nothing, exchange=nothing)::F where {F<:Float64} + +Run IHACRESBilinearNodeGW for one timestep. + +# Arguments +- `node` : IHACRESBilinearNodeGW instance +- `rain` : Rainfall (mm/day) +- `evap` : Evaporation or temperature +- `ts` : Current timestep index +- `inflow` : Inflow from upstream nodes (ML/day), default 0.0 +- `extraction` : Groundwater extraction (ML/day), default 0.0 +- `exchange` : External GW flux (ML/day), default 0.0 + +# Returns +Total outflow (ML/day) = quickflow + baseflow + inflow - extraction +""" +function run_timestep!( + node::IHACRESBilinearNodeGW, + rain::F, + evap::F, + ts::Int64; + inflow=nothing, + extraction=nothing, + exchange=nothing +)::F where {F<:Float64} + # Get current states + current_cmd::F = node.storage[ts] + quick_store::F = node.quick_store[ts] + gw_store::F = node.gw_storage[ts] + + # Parse optional inputs + in_flow::F = timestep_value(ts, node.name, "inflow", inflow) + gw_ext::F = timestep_value(ts, node.name, "extraction", extraction) + gw_exchange::F = timestep_value(ts, node.name, "exchange", exchange) + + # 1. Calculate interim CMD and effective rainfall + (mf, e_rainfall, recharge) = calc_ft_interim_cmd( + current_cmd, + rain, + node.d.val::F, + node.d2.val::F, + node.alpha.val::F + ) + + # 2. Calculate ET + et::F = calc_ET_from_E( + node.e.val::F, + evap, + mf, + node.f.val::F, + node.d.val::F + ) + + # 3. Update CMD + cmd::F = calc_cmd( + current_cmd, + rain, + et, + e_rainfall, + recharge + ) + + # 4. Calculate quickflow (uses 1-vs fraction of effective rainfall) + (nq_store, quickflow) = calc_ft_quick_flow( + quick_store, + e_rainfall, + node.area, + node.a.val::F, + node.vs.val::F + ) + + # 5. Calculate GW storage and baseflow (uses vs fraction of effective rainfall) + (ngw_store, baseflow) = routing( + gw_store, + e_rainfall, + node.area, + node.tau_s.val::F, + node.vs.val::F, + gw_ext, + node.L.val::F, + gw_exchange + ) + + # 6. Calculate total outflow + # Total = local quickflow + local baseflow + upstream inflow - local extraction + outflow = quickflow + baseflow + in_flow + + # 7. Update state + node.storage[ts+1] = cmd + node.quick_store[ts+1] = nq_store + node.gw_storage[ts+1] = ngw_store + + node.inflow[ts] = in_flow + node.quickflow[ts] = quickflow + node.baseflow[ts] = baseflow + node.outflow[ts] = outflow + node.effective_rainfall[ts] = e_rainfall + node.et[ts] = et + + return outflow +end + + +""" + param_info(node::IHACRESBilinearNodeGW)::Tuple + +Extract parameter names, values, and bounds for calibration. + +Returns tuple of (param_names, values, bounds). +""" +function param_info(node::IHACRESBilinearNodeGW)::Tuple + tmp = Model(node) + values = collect(tmp[:val]) + bounds = collect(tmp[:bounds]) + param_names = collect(tmp[:fieldname]) + + return param_names, values, bounds +end + + +""" + update_params!(node::IHACRESBilinearNodeGW, d, d2, e, f, a, alpha, tau_s, vs, L, specific_yield) + +Update node parameters (typically after calibration). +""" +function update_params!( + node::IHACRESBilinearNodeGW, + d::Float64, d2::Float64, e::Float64, f::Float64, + a::Float64, alpha::Float64, + tau_s::Float64, vs::Float64, L::Float64, specific_yield::Float64 +)::Nothing + node.d = Param(d, bounds=node.d.bounds) + node.d2 = Param(d2, bounds=node.d2.bounds) + node.e = Param(e, bounds=node.e.bounds) + node.f = Param(f, bounds=node.f.bounds) + node.a = Param(a, bounds=node.a.bounds) + node.alpha = Param(alpha, bounds=node.alpha.bounds) + node.tau_s = Param(tau_s, bounds=node.tau_s.bounds) + node.vs = Param(vs, bounds=node.vs.bounds) + node.L = Param(L, bounds=node.L.bounds) + node.specific_yield = Param(specific_yield, bounds=node.specific_yield.bounds) + + return nothing +end + + +""" + extract_spec!(node::IHACRESBilinearNodeGW, spec::AbstractDict)::Nothing + +Extract node-specific details to specification dictionary. +""" +function extract_spec!(node::IHACRESBilinearNodeGW, spec::AbstractDict)::Nothing + spec["initial_storage"] = node.storage[1] + spec["initial_gw_storage"] = node.gw_storage[1] + spec["aquifer_area"] = node.aquifer_area + spec["gauge_elevation"] = node.gauge_elevation + spec["bore_ground_elevation"] = node.bore_ground_elevation + + return nothing +end + + +""" + get_simulated_depth(node::IHACRESBilinearNodeGW, ts::Int64)::Float64 + +Get simulated depth below ground surface for comparison with bore observations. + +Returns depth in meters below ground surface at bore location. +Positive values indicate water table is below ground surface (normal condition). +""" +function get_simulated_depth(node::IHACRESBilinearNodeGW, ts::Int64)::Float64 + return convert_storage_to_depth( + node.gw_storage[ts], + node.aquifer_area, # Use aquifer_area, not catchment area + node.specific_yield.val, + node.gauge_elevation, + node.bore_ground_elevation + ) +end + + +""" + get_water_table_elevation(node::IHACRESBilinearNodeGW, ts::Int64)::Float64 + +Get water table elevation in mAHD at timestep ts. + +Returns absolute elevation of water table in Australian Height Datum (mAHD). +This is the elevation of the water table surface, useful for regional mapping +and comparing water levels between different locations. + +# Example +```julia +wt_elevation = get_water_table_elevation(node, 100) # Get elevation at day 100 +println("Water table at ", wt_elevation, " mAHD") +``` +""" +function get_water_table_elevation(node::IHACRESBilinearNodeGW, ts::Int64)::Float64 + area_m2 = node.aquifer_area * 1e6 # Use aquifer_area + water_table_height = (node.gw_storage[ts] * 1000.0) / (area_m2 * node.specific_yield.val) + return node.gauge_elevation + water_table_height +end + + +""" + get_water_table_height_above_streambed(node::IHACRESBilinearNodeGW, ts::Int64)::Float64 + +Get height of water table above stream bed in meters at timestep ts. + +Returns height above gauge elevation (the reference where storage=0). +Positive values indicate water table is above stream bed (normal for gaining stream). +Negative values indicate water table is below stream bed (losing stream or disconnected). + +# Example +```julia +height = get_water_table_height_above_streambed(node, 100) +if height > 0 + println("Water table is ", height, " m above stream bed - baseflow likely") +else + println("Water table is below stream bed - no baseflow") +end +``` +""" +function get_water_table_height_above_streambed(node::IHACRESBilinearNodeGW, ts::Int64)::Float64 + area_m2 = node.aquifer_area * 1e6 # Use aquifer_area + return (node.gw_storage[ts] * 1000.0) / (area_m2 * node.specific_yield.val) +end diff --git a/src/Nodes/IHACRES/components/gw_flow.jl b/src/Nodes/IHACRES/components/gw_flow.jl new file mode 100644 index 00000000..46ba2927 --- /dev/null +++ b/src/Nodes/IHACRES/components/gw_flow.jl @@ -0,0 +1,248 @@ +""" + routing( + prev_gw_storage::Float64, + recharge::Float64, + area::Float64, + tau_s::Float64, + vs::Float64, + extraction::Float64, + loss::Float64, + gw_exchange::Float64 + )::Tuple{Float64,Float64} + +Calculate groundwater storage and baseflow using IHACRES_GW formulation. + +Based on Ivkovic et al. (2005) groundwater module extension to IHACRES. + +# Arguments +- `prev_gw_storage` : Groundwater storage at t-1 (ML) +- `recharge` : Recharge from effective rainfall (mm/day) +- `area` : Catchment area (km²) +- `tau_s` : Slow flow time constant (days), controls recession rate +- `vs` : Recharge fraction (0-1), proportion of effective rainfall to GW +- `extraction` : Groundwater extraction/pumping (ML/day) +- `loss` : Constant loss term L (ML/day), positive = loss, negative = gain +- `gw_exchange` : External GW flux (ML/day), positive = inflow, negative = outflow + +# Returns +Tuple of (gw_storage, baseflow) both in ML/day + +# References +- Ivkovic et al. (2005a) - IHACRES_GW model description +- Croke & Jakeman (2004) - IHACRES-CMD module + +# Notes +- When gw_storage ≤ 0, baseflow = 0 (water table below stream bed) +- The coefficient 'a' is derived from tau_s: a = 1/τ_s (exponential recession coefficient) +- Recharge is converted from mm/day to ML/day by: recharge_ML = vs * recharge * area +- Storage and baseflow are partitioned as: G[t] = G_tmp/(1+a), Q[t] = a*G_tmp/(1+a) +""" +function routing( + prev_gw_storage::Float64, + recharge::Float64, + area::Float64, + tau_s::Float64, + vs::Float64, + extraction::Float64, + loss::Float64, + gw_exchange::Float64 +)::Tuple{Float64,Float64} + # Calculate baseflow coefficient 'a' from time constant tau_s + # For exponential recession: a = 1/τ_s + # This gives Q = a*G which produces exponential decay + a = 1.0 / tau_s + + # Convert recharge from mm/day to ML/day + # vs is the fraction that goes to groundwater (remainder goes to quick flow) + recharge_ML = vs * recharge * area + + # Update storage (before baseflow calculation) + # Storage increases with recharge and gw_exchange (inflow) + # Storage decreases with extraction and loss + gw_tmp = prev_gw_storage + recharge_ML - extraction - loss + gw_exchange + + # Calculate baseflow and final storage + # Baseflow only occurs if storage is positive (water table above stream bed) + if gw_tmp > 0.0 + # Split storage into baseflow (discharge) and remaining storage + # From IHACRES_GW: Q_s = a * G_t and G_t = G_tmp / (1 + a) + baseflow = a * gw_tmp / (1.0 + a) + gw_storage = gw_tmp / (1.0 + a) + else + # Water table below stream bed - no baseflow + baseflow = 0.0 + gw_storage = max(0.0, gw_tmp) # Prevent negative storage + end + + # Ensure baseflow is non-negative + baseflow = max(0.0, baseflow) + + return (gw_storage, baseflow) +end + + +""" + calc_ft_quick_flow( + prev_quick::Float64, + e_rain::Float64, + area::Float64, + a::Float64, + vs::Float64 + )::Tuple{Float64,Float64} + +Calculate quickflow store and discharge. + +Extracted from calc_ft_flows to allow separate quick and slow flow calculation +in IHACRES_GW where slow flow is replaced by groundwater module. + +# Arguments +- `prev_quick` : Previous quickflow storage (ML) +- `e_rain` : Effective rainfall (mm/day) +- `area` : Catchment area (km²) +- `a` : Quickflow storage coefficient (1/τ_q) +- `vs` : Recharge fraction (proportion going to GW, remainder goes to quick) + +# Returns +Tuple of (quick_store, quickflow) in ML/day + +# Notes +- Only (1 - vs) fraction of effective rainfall goes to quickflow +- vs fraction goes to groundwater recharge instead +""" +function calc_ft_quick_flow( + prev_quick::Float64, + e_rain::Float64, + area::Float64, + a::Float64, + vs::Float64 +)::Tuple{Float64,Float64} + # Partition effective rainfall between quickflow and GW recharge + # (1 - vs) goes to quick flow, vs goes to GW recharge + quick_rain = (1.0 - vs) * e_rain + + # Calculate quickflow store and discharge + tmp_calc = max(0.0, prev_quick + (quick_rain * area)) + + if tmp_calc > 0.0 + alpha = exp(-a) + beta = (1.0 - alpha) * tmp_calc + quick_store = alpha * tmp_calc + quickflow = beta + else + quick_store = tmp_calc + quickflow = 0.0 + end + + @assert quickflow >= 0.0 "Quickflow cannot be negative: $quickflow" + + return (quick_store, quickflow) +end + + +""" + convert_storage_to_depth( + gw_storage_ML::Float64, + area_km2::Float64, + specific_yield::Float64, + stream_bed_elevation_mAHD::Float64, + bore_ground_elevation_mAHD::Float64 + )::Float64 + +Convert groundwater storage (ML) to depth below ground surface (m). + +Used for comparing simulated water table with bore observations of "Depth to water (m)". + +# Arguments +- `gw_storage_ML` : Groundwater storage volume (ML) +- `area_km2` : Catchment area (km²) +- `specific_yield` : Aquifer specific yield (dimensionless, typically 0.1-0.3) +- `stream_bed_elevation_mAHD` : Stream bed elevation (mAHD), reference where storage=0 +- `bore_ground_elevation_mAHD` : Ground surface elevation at bore (mAHD) + +# Returns +Depth below ground surface (m), positive values indicate water table below ground + +# Calculation Steps +1. Convert storage to water table height above stream bed: + height = storage (ML) * 1000 / (area_m² * specific_yield) +2. Calculate water table elevation: + elevation = stream_bed + height +3. Calculate depth below surface: + depth = ground_surface - water_table_elevation + +# Notes +- When storage = 0: water table is at stream bed level +- Positive depth = water table below ground (normal condition) +- Negative depth = water table above ground (flooding/artesian) +""" +function convert_storage_to_depth( + gw_storage_ML::Float64, + area_km2::Float64, + specific_yield::Float64, + stream_bed_elevation_mAHD::Float64, + bore_ground_elevation_mAHD::Float64 +)::Float64 + # Convert area from km² to m² + area_m2 = area_km2 * 1e6 + + # Convert storage (ML) to water table height (m) above stream bed + # Storage in ML → m³ (×1000) → divide by (area × specific_yield) + water_table_height_m = (gw_storage_ML * 1000.0) / (area_m2 * specific_yield) + + # Calculate water table elevation in mAHD + # When storage = 0, water table is at stream bed elevation + water_table_elevation_mAHD = stream_bed_elevation_mAHD + water_table_height_m + + # Calculate depth below ground surface (positive = below ground) + depth_below_surface_m = bore_ground_elevation_mAHD - water_table_elevation_mAHD + + return depth_below_surface_m +end + + +""" + convert_depth_to_storage( + depth_below_surface_m::Float64, + area_km2::Float64, + specific_yield::Float64, + stream_bed_elevation_mAHD::Float64, + bore_ground_elevation_mAHD::Float64 + )::Float64 + +Convert depth below ground surface (m) to groundwater storage (ML). + +Inverse of convert_storage_to_depth. Useful for initializing storage from observed bore depths. + +# Arguments +- `depth_below_surface_m` : Observed depth to water below ground (m) +- `area_km2` : Catchment area (km²) +- `specific_yield` : Aquifer specific yield (dimensionless) +- `stream_bed_elevation_mAHD` : Stream bed elevation (mAHD) +- `bore_ground_elevation_mAHD` : Ground surface elevation at bore (mAHD) + +# Returns +Groundwater storage (ML) + +# Notes +- Used to back-calculate initial storage from first bore observation +- Can return negative storage if water table is below stream bed +""" +function convert_depth_to_storage( + depth_below_surface_m::Float64, + area_km2::Float64, + specific_yield::Float64, + stream_bed_elevation_mAHD::Float64, + bore_ground_elevation_mAHD::Float64 +)::Float64 + # Calculate water table elevation from depth + water_table_elevation_mAHD = bore_ground_elevation_mAHD - depth_below_surface_m + + # Calculate height above stream bed + water_table_height_m = water_table_elevation_mAHD - stream_bed_elevation_mAHD + + # Convert to storage + area_m2 = area_km2 * 1e6 + gw_storage_ML = (water_table_height_m * area_m2 * specific_yield) / 1000.0 + + return gw_storage_ML +end diff --git a/src/Nodes/Node.jl b/src/Nodes/Node.jl index f6138a7a..88c60785 100644 --- a/src/Nodes/Node.jl +++ b/src/Nodes/Node.jl @@ -162,11 +162,10 @@ function Base.show(io::IO, n::NetworkNode) param_names = n_model[:fieldname] x0 = n_model[:val] bounds = n_model[:bounds] - descs = n_model[:desc] lb, ub = zip(bounds...) - details = hcat([param_names...], [x0...], [lb...], [ub...], [descs...]) + details = hcat([param_names...], [x0...], [lb...], [ub...]) - pretty_table(io, details, header=["Parameter", "Value", "Lower Bound", "Upper Bound", "Description"]) + pretty_table(io, details, header=["Parameter", "Value", "Lower Bound", "Upper Bound"]) print("\n") end diff --git a/src/Streamfall.jl b/src/Streamfall.jl index 0c30579e..cd6d1085 100644 --- a/src/Streamfall.jl +++ b/src/Streamfall.jl @@ -12,6 +12,7 @@ include("calibration.jl") include("Analysis/Analysis.jl") include("Nodes/IHACRES/IHACRESNode.jl") +include("Nodes/IHACRES/IHACRESBilinearNodeGW.jl") # include("Nodes/IHACRES/IHACRESExpuhNode.jl") include("Nodes/GR4J/GR4JNode.jl") include("Nodes/HyMod/HyModNode.jl") @@ -278,7 +279,7 @@ end # Nodes export NetworkNode, GenericNode, GenericDirectNode -export IHACRES, IHACRESNode, IHACRESBilinearNode, DamNode, Climate +export IHACRES, IHACRESNode, IHACRESBilinearNode, IHACRESBilinearNodeGW, DamNode, Climate export create_node, GR4JNode, HyModNode, SimpleHyModNode, SIMHYDNode export EnsembleNode, WeightedEnsembleNode, GREnsembleNode @@ -310,4 +311,10 @@ export Analysis export apply_bias_correction +# IHACRES_GW functions (for testing and utilities) +export convert_storage_to_depth, convert_depth_to_storage +export calc_ft_quick_flow, routing +export get_simulated_depth +export get_water_table_elevation, get_water_table_height_above_streambed + end # end module diff --git a/test/data/campaspe/campaspe_data_prep.jl b/test/data/campaspe/campaspe_data_prep.jl index 481142ce..15ad06cb 100644 --- a/test/data/campaspe/campaspe_data_prep.jl +++ b/test/data/campaspe/campaspe_data_prep.jl @@ -16,6 +16,8 @@ climate_data = CSV.read( comment="#" ) +sn = load_network("Example Network", "campaspe_network.yml") + # Load in observation data for each gauge in network outflow_files = readdir(glob"gauges/*_outflow*.csv") all_flow_data = CSV.read.(outflow_files, DataFrame; comment="#") @@ -27,11 +29,11 @@ for name in node_names(sn) suffix = "_outflow_[ML]" if name != "406000" - ordered_flow_data[name] = extract_flow(all_flow_data[matching_file_pos], name, suffix) + ordered_flow_data[name] = extract_flow(all_flow_data[matching_file_pos], name; suffix=suffix) continue end - out = extract_flow(all_flow_data[matching_file_pos], name, suffix) + out = extract_flow(all_flow_data[matching_file_pos], name;suffix=suffix) ext = all_flow_data[matching_file_pos][:, ["Date", "406000_extractions_[ML]"]] df = DataFrame(Date=out.Date) diff --git a/test/runtests.jl b/test/runtests.jl index ebaa6dd9..c000c36c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -174,8 +174,9 @@ end # @test flow_results[2] == slow_store # end -include("test_metrics.jl") -include("test_data_op.jl") -include("test_nodes.jl") -include("test_networks.jl") -include("test_calibration.jl") +#include("test_metrics.jl") +#include("test_data_op.jl") +#include("test_nodes.jl") +#include("test_networks.jl") +#include("test_calibration.jl") +include("test_ihacres_gw.jl") diff --git a/test/test_ihacres_gw.jl b/test/test_ihacres_gw.jl new file mode 100644 index 00000000..5815d554 --- /dev/null +++ b/test/test_ihacres_gw.jl @@ -0,0 +1,461 @@ +using Test +using Streamfall + + +@testset "IHACRES_GW Core Functions" begin + + @testset "routing" begin + # Test basic GW storage and baseflow calculation + + # Parameters + prev_storage = 1000.0 # ML + recharge = 5.0 # mm/day + area = 254.0 # km² + tau_s = 20.0 # days + vs = 0.5 # fraction to GW + extraction = 10.0 # ML/day + loss = 0.0 # ML/day + gw_exchange = 0.0 # ML/day + + (new_storage, baseflow) = routing( + prev_storage, recharge, area, tau_s, vs, extraction, loss, gw_exchange + ) + + # Check that storage and baseflow are non-negative + @test new_storage >= 0.0 + @test baseflow >= 0.0 + + # Check water balance (approximately) + recharge_ML = vs * recharge * area + # Storage change should approximately equal: inflows - outflows + # new_storage ≈ prev_storage + recharge_ML - baseflow - extraction + # (approximately, because baseflow depends on storage) + storage_change = new_storage - prev_storage + net_flux = recharge_ML - baseflow - extraction + @test abs(storage_change - net_flux) < 1.0 # Within 1 ML tolerance + + # Test zero storage case (no baseflow) + (zero_storage, zero_baseflow) = routing( + 0.0, 0.0, area, tau_s, vs, 0.0, 0.0, 0.0 + ) + @test zero_storage == 0.0 + @test zero_baseflow == 0.0 + + # Test negative storage input (should be clamped to zero) + (neg_storage, neg_baseflow) = routing( + -100.0, 0.0, area, tau_s, vs, 0.0, 0.0, 0.0 + ) + @test neg_storage == 0.0 # Storage clamped to zero (cannot be negative) + @test neg_baseflow == 0.0 # No baseflow when storage is zero + end + + + @testset "calc_ft_quick_flow" begin + # Test quickflow calculation + + prev_quick = 100.0 # ML + e_rain = 10.0 # mm/day + area = 254.0 # km² + a = 0.9 # quickflow coefficient + vs = 0.5 # fraction to GW (so 0.5 to quick) + + (new_quick, quickflow) = calc_ft_quick_flow(prev_quick, e_rain, area, a, vs) + + # Check non-negative + @test new_quick >= 0.0 + @test quickflow >= 0.0 + + # Check that quickflow uses (1-vs) fraction + expected_quick_rain = (1.0 - vs) * e_rain + # new_quick should be based on expected_quick_rain + @test new_quick > 0.0 # Should have some storage + + # Test zero case + (zero_quick, zero_qflow) = calc_ft_quick_flow(0.0, 0.0, area, a, vs) + @test zero_quick == 0.0 + @test zero_qflow == 0.0 + end + + + @testset "convert_storage_to_depth" begin + # Test storage to depth conversion + + storage = 5000.0 # ML + area = 254.0 # km² + sy = 0.15 # specific yield + stream_bed = 120.0 # mAHD + bore_ground = 123.0 # mAHD + + depth = convert_storage_to_depth(storage, area, sy, stream_bed, bore_ground) + + # Depth should be positive (water table below ground) + @test depth > 0.0 + + # Manual calculation check + area_m2 = area * 1e6 + wt_height = (storage * 1000.0) / (area_m2 * sy) + wt_elev = stream_bed + wt_height + expected_depth = bore_ground - wt_elev + @test abs(depth - expected_depth) < 1e-6 + + # Test zero storage (water table at stream bed) + depth_zero = convert_storage_to_depth(0.0, area, sy, stream_bed, bore_ground) + @test abs(depth_zero - (bore_ground - stream_bed)) < 1e-6 + + # Test negative depth (water table above ground - unlikely but possible) + large_storage = 50000.0 # ML + depth_neg = convert_storage_to_depth(large_storage, area, sy, stream_bed, bore_ground) + # Could be negative if storage is very large + end + + + @testset "convert_depth_to_storage" begin + # Test depth to storage conversion (inverse operation) + + depth = 10.0 # m below ground + area = 254.0 # km² + sy = 0.15 + stream_bed = 120.0 # mAHD + bore_ground = 123.0 # mAHD + + storage = convert_depth_to_storage(depth, area, sy, stream_bed, bore_ground) + + # Storage could be negative if depth is large (water table below stream) + # But should be calculable + @test !isnan(storage) + + # Test round-trip conversion + storage_test = 5000.0 # ML + depth_from_storage = convert_storage_to_depth( + storage_test, area, sy, stream_bed, bore_ground + ) + storage_back = convert_depth_to_storage( + depth_from_storage, area, sy, stream_bed, bore_ground + ) + @test abs(storage_back - storage_test) < 1e-3 # Within 1e-3 ML + end + +end + + +@testset "IHACRESBilinearNodeGW Node" begin + + @testset "Node construction from dict" begin + # Test creating node from specification dictionary + + spec = Dict( + "area" => 254.0, + "aquifer_area" => 254.0, + "initial_storage" => 100.0, + "initial_gw_storage" => 1000.0, + "gauge_elevation" => 120.0, + "bore_ground_elevation" => 123.0, + "parameters" => Dict( + "d" => 200.0, + "d2" => 2.0, + "e" => 1.0, + "f" => 0.8, + "a" => 0.9, + "alpha" => 0.1, + "tau_s" => 20.0, + "vs" => 0.5, + "L" => 0.0, + "specific_yield" => 0.15 + ) + ) + + node = IHACRESBilinearNodeGW("test_node", spec) + + @test node.name == "test_node" + @test node.area == 254.0 + @test node.storage[1] == 100.0 + @test node.gw_storage[1] == 1000.0 + @test node.gauge_elevation == 120.0 + @test node.bore_ground_elevation == 123.0 + @test node.tau_s.val == 20.0 + @test node.vs.val == 0.5 + end + + + @testset "Node construction direct" begin + # Test creating node with direct parameters + + node = IHACRESBilinearNodeGW( + "test_direct", + 254.0, # area + 200.0, 2.0, 1.0, 0.8, # d, d2, e, f + 0.9, 0.1, # a, alpha + 20.0, 0.5, 0.0, 0.15, # tau_s, vs, L, sy + 120.0, 123.0, # gauge_elev, bore_elev + 100.0, 0.0, 1000.0 # cmd, quick, gw stores + ) + + @test node.name == "test_direct" + @test node.area == 254.0 + @test node.tau_s.val == 20.0 + end + + + @testset "prep_state!" begin + # Test state preparation + + node = IHACRESBilinearNodeGW( + "test_prep", + 254.0, + 200.0, 2.0, 1.0, 0.8, + 0.9, 0.1, + 20.0, 0.5, 0.0, 0.15, + 120.0, 123.0, + 100.0, 0.0, 1000.0 + ) + + timesteps = 365 + prep_state!(node, timesteps) + + @test length(node.storage) == timesteps + 1 + @test length(node.gw_storage) == timesteps + 1 + @test length(node.quick_store) == timesteps + 1 + @test length(node.outflow) == timesteps + @test length(node.baseflow) == timesteps + @test length(node.quickflow) == timesteps + end + + + @testset "reset!" begin + # Test reset functionality + + node = IHACRESBilinearNodeGW( + "test_reset", + 254.0, + 200.0, 2.0, 1.0, 0.8, + 0.9, 0.1, + 20.0, 0.5, 0.0, 0.15, + 120.0, 123.0, + 100.0, 0.0, 1000.0 + ) + + prep_state!(node, 10) + + # Run a few timesteps to modify state + node.storage[2] = 90.0 + node.gw_storage[2] = 1100.0 + + reset!(node) + + @test length(node.storage) == 1 + @test length(node.gw_storage) == 1 + @test node.storage[1] == 100.0 # Back to initial + @test node.gw_storage[1] == 1000.0 # Back to initial + end + + + @testset "run_timestep!" begin + # Test running a single timestep + + node = IHACRESBilinearNodeGW( + "test_run", + 254.0, + 200.0, 2.0, 1.0, 0.8, + 0.9, 0.1, + 20.0, 0.5, 0.0, 0.15, + 120.0, 123.0, + 100.0, 0.0, 1000.0 + ) + + prep_state!(node, 1) + + rain = 10.0 # mm + evap = 5.0 # mm + ts = 1 + + outflow = run_timestep!(node, rain, evap, ts) + + # Check that outflow was calculated + @test outflow >= 0.0 + @test node.outflow[1] == outflow + + # Check that states were updated + @test node.storage[2] != node.storage[1] + @test node.quickflow[1] >= 0.0 + @test node.baseflow[1] >= 0.0 + + # Check that quickflow + baseflow contributes to outflow + @test node.outflow[1] >= node.quickflow[1] + node.baseflow[1] + end + + + @testset "get_simulated_depth" begin + # Test getting simulated depth for calibration + + node = IHACRESBilinearNodeGW( + "test_depth", + 254.0, + 200.0, 2.0, 1.0, 0.8, + 0.9, 0.1, + 20.0, 0.5, 0.0, 0.15, + 120.0, 123.0, + 100.0, 0.0, 5000.0 # 5000 ML GW storage + ) + + depth = get_simulated_depth(node, 1) + + @test depth > 0.0 # Should be positive (below ground) + @test !isnan(depth) + + # Check consistency with convert_storage_to_depth + expected_depth = convert_storage_to_depth( + node.gw_storage[1], + node.area, + node.specific_yield.val, + node.gauge_elevation, + node.bore_ground_elevation + ) + @test abs(depth - expected_depth) < 1e-6 + end + +end + + +@testset "Parameter extraction" begin + + node = IHACRESBilinearNodeGW( + "test_params", + 254.0, + 200.0, 2.0, 1.0, 0.8, + 0.9, 0.1, + 20.0, 0.5, 0.0, 0.15, + 120.0, 123.0, + 100.0, 0.0, 1000.0 + ) + + param_names, values, bounds = param_info(node) + + @test length(param_names) == length(values) + @test length(param_names) == length(bounds) + @test length(param_names) == 10 # 10 calibratable parameters + + # Check that key parameters are present + @test :tau_s in param_names + @test :vs in param_names + @test :L in param_names + @test :specific_yield in param_names + +end + + +@testset "Water Balance Partitioning" begin + # Test that effective rainfall is correctly partitioned between quickflow and GW + # This is critical for ensuring the model conserves water properly + + @testset "Effective rainfall partitioning" begin + node = IHACRESBilinearNodeGW( + "test_balance", + 100.0, # area + 50.0, 2.0, 1.0, 0.8, # d, d2, e, f (low d to generate effective rainfall) + 0.5, 0.5, # a, alpha + 20.0, 0.6, 0.0, 0.15, # tau_s, vs=0.6, L=0, sy + 120.0, 123.0, + 10.0, 0.0, 1000.0 # low initial CMD to generate effective rainfall + ) + + prep_state!(node, 10) + + # Run a few timesteps with significant rainfall + rain_events = [50.0, 30.0, 40.0, 20.0, 10.0] + evap = 5.0 + + for (i, rain) in enumerate(rain_events) + run_timestep!(node, rain, evap, i) + end + + # Check water balance for days with effective rainfall + vs = node.vs.val + area = node.area + + for t in 1:length(rain_events) + eff_rain = node.effective_rainfall[t] + + if eff_rain > 0.01 # Only check days with significant effective rainfall + # Convert effective rainfall to ML + eff_rain_ML = eff_rain * area + + # Check quickflow + # Note: quickflow includes previous storage, so we need to account for that + # Expected addition to quick from this timestep: (1-vs) * eff_rain * area + expected_quick_rain_ML = (1.0 - vs) * eff_rain * area + + # Check GW recharge + # Expected addition to GW from this timestep: vs * eff_rain * area + # But the code uses: vs * recharge * area, where recharge != eff_rain + expected_gw_rain_ML = vs * eff_rain * area + + # Calculate actual GW inflow (storage change + baseflow + L + extraction) + storage_change = node.gw_storage[t+1] - node.gw_storage[t] + baseflow = node.baseflow[t] + L = node.L.val + actual_gw_inflow = storage_change + baseflow + L # No extraction in this test + + # Print diagnostic info + println("\nTimestep $t:") + println(" Effective rainfall: $(round(eff_rain, digits=2)) mm = $(round(eff_rain_ML, digits=2)) ML") + println(" Expected to quick ($(round(1-vs, digits=2))*eff_rain): $(round(expected_quick_rain_ML, digits=2)) ML") + println(" Expected to GW ($(round(vs, digits=2))*eff_rain): $(round(expected_gw_rain_ML, digits=2)) ML") + println(" Actual to GW (storage change + baseflow + L): $(round(actual_gw_inflow, digits=2)) ML") + println(" Difference (actual - expected): $(round(actual_gw_inflow - expected_gw_rain_ML, digits=2)) ML") + println(" Ratio (actual/expected): $(round(actual_gw_inflow / expected_gw_rain_ML, digits=3))") + + # TEST: The partitioning should satisfy vs * eff_rain going to GW + # Allow 5% tolerance for numerical issues + @test isapprox(actual_gw_inflow, expected_gw_rain_ML, rtol=0.05) + end + end + end + + @testset "Total water conservation" begin + # Test that water is conserved over entire simulation + node = IHACRESBilinearNodeGW( + "test_conservation", + 100.0, + 50.0, 2.0, 1.0, 0.8, + 0.5, 0.5, + 20.0, 0.5, 0.0, 0.15, # vs=0.5 for equal split + 120.0, 123.0, + 10.0, 0.0, 1000.0 + ) + + prep_state!(node, 20) + + # Run simulation + for t in 1:20 + rain = 20.0 + 10.0 * sin(t / 3.0) # Variable rainfall + evap = 5.0 + run_timestep!(node, rain, evap, t) + end + + # Sum all water fluxes + total_eff_rain = sum(node.effective_rainfall) * node.area + total_quickflow = sum(node.quickflow) + total_baseflow = sum(node.baseflow) + + # Storage changes + gw_storage_change = node.gw_storage[end] - node.gw_storage[1] + quick_storage_change = node.quick_store[end] - node.quick_store[1] + + # Total water out + storage change should equal total effective rainfall + total_out = total_quickflow + total_baseflow + total_storage_change = gw_storage_change + quick_storage_change + + println("\nWater Balance Summary:") + println(" Total effective rainfall: $(round(total_eff_rain, digits=1)) ML") + println(" Total quickflow: $(round(total_quickflow, digits=1)) ML") + println(" Total baseflow: $(round(total_baseflow, digits=1)) ML") + println(" GW storage change: $(round(gw_storage_change, digits=1)) ML") + println(" Quick storage change: $(round(quick_storage_change, digits=1)) ML") + println(" Total out + storage change: $(round(total_out + total_storage_change, digits=1)) ML") + println(" Difference: $(round(total_eff_rain - (total_out + total_storage_change), digits=1)) ML") + println(" Percentage: $(round(100*(total_out + total_storage_change)/total_eff_rain, digits=1))%") + + # Water balance check (allow 1% error) + @test isapprox(total_eff_rain, total_out + total_storage_change, rtol=0.01) + end +end