diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 5592cc7c..f0a02876 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -13,7 +13,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.10' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - '1.11' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
           - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
           - 'nightly'
         os:
diff --git a/Project.toml b/Project.toml
index f59ec6bf..6555a2d5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -52,6 +52,7 @@ CodecZlib = "0.7"
 Compose = "0.9"
 DataFrames = "1.6"
 DataStructures = "0.18.10"
+Dates = "1.11.0"
 Distributions = "0.25.100"
 GR = "0.72.9, 0.73"
 GR_jll = "0.73"
diff --git a/README.md b/README.md
index 82819cde..54452715 100644
--- a/README.md
+++ b/README.md
@@ -9,7 +9,7 @@
 Streamfall leverages the Julia language and ecosystem to provide:
 - Quick heterogeneous modelling of a stream network
 - Use of different rainfall-runoff models and their ensembles in tandem
-- Modelling and assessment of interacting systems
+- Support for modelling and assessment of interacting systems
 - A wide range of performance metrics
 
 This package includes implementations of the following:
@@ -20,12 +20,21 @@ This package includes implementations of the following:
 
 Performance is expected to be similar to implementations in C and Fortran.
 
+Naive timings (using `@time`) for an example dataset spanning 1963-07-05 - 2014-12-31 (18808 days, approximately 51.5 years)
+
+- SimpleHyModNode \
+  0.016502 seconds (469.25 k allocations: 12.902 MiB)
+- GR4JNode \
+  0.015274 seconds (224.75 k allocations: 5.584 MiB)
+- SYMHYDNode \
+  0.039540 seconds (638.01 k allocations: 15.190 MiB, 46.99% gc time)
+- IHACRESBilinearNode \
+  0.021734 seconds (675.63 k allocations: 17.773 MiB)
+
 The IHACRES rainfall-runoff model was previously implemented with [ihacres_nim](https://github.com/ConnectedSystems/ihacres_nim) but has since been ported to pure Julia.
 
 [Graphs](https://github.com/JuliaGraphs/Graphs.jl) and [MetaGraphs](https://github.com/JuliaGraphs/MetaGraphs.jl) are used underneath for network traversal/analysis.
 
-> [NOTE] Streamfall is currently in its early stages and under active development. Although it is fairly usable for small networks and assessing/analyzing single sub-catchments, things may change drastically and unexpectedly.
-
 ## Installation
 
 Streamfall is now registered! The latest release version can be installed with:
diff --git a/docs/src/API/nodes/IHACRES.md b/docs/src/API/nodes/IHACRES.md
index 681f24ba..4a18bfc4 100644
--- a/docs/src/API/nodes/IHACRES.md
+++ b/docs/src/API/nodes/IHACRES.md
@@ -5,11 +5,3 @@ Modules = [Streamfall]
 Order   = [:function, :type]
 Pages   = ["Nodes/IHACRES/IHACRESNode.jl"]
 ```
-
-### IHACRES - Expuh
-
-```@autodocs
-Modules = [Streamfall]
-Order   = [:function, :type]
-Pages   = ["IHACRESExpuhNode.jl"]
-```
diff --git a/examples/node_comparison.jl b/examples/node_comparison.jl
index e2e24595..5cac290b 100644
--- a/examples/node_comparison.jl
+++ b/examples/node_comparison.jl
@@ -1,4 +1,6 @@
-using Distributed, Plots, StatsPlots
+using Distributed
+using Plots
+
 
 N = 4
 if nworkers() < N
@@ -18,7 +20,7 @@ end
         comment="#",
         dateformat=date_format) |> DataFrame
 
-    hist_streamflow = obs_data[:, ["Date", "410730_Q"]]
+    hist_streamflow = extract_flow(obs_data, "410730")
     climate_data = obs_data[:, ["Date", "410730_P", "410730_PET"]]
     climate = Climate(climate_data, "_P", "_PET")
 
@@ -26,7 +28,7 @@ end
 
     # Create objective function to minimize (here we use Normalized KGE')
     func = (obs, sim) -> 1.0 - Streamfall.NmKGE(obs[burn_in:end], sim[burn_in:end])
-    opt_func = (node) -> calibrate!(node, climate, hist_streamflow[:, "410730_Q"]; metric=func, MaxTime=900)
+    opt_func = (node) -> calibrate!(node, climate, hist_streamflow[:, "410730"], func; MaxTime=30)
 end
 
 
@@ -34,7 +36,7 @@ end
 hymod_node = create_node(SimpleHyModNode, "410730", 129.2)
 gr4j_node = create_node(GR4JNode, "410730", 129.2)
 symhyd_node = create_node(SYMHYDNode, "410730", 129.2)
-ihacres_node = create_node(BilinearNode, "410730", 129.2)
+ihacres_node = create_node(IHACRESBilinearNode, "410730", 129.2)
 
 
 # Calibrate each node separately using multiprocessing
@@ -44,14 +46,14 @@ result = pmap(opt_func, node_list)
 
 
 # Create comparison plot
-Qo = hist_streamflow[:, "410730_Q"]
-Qo_burn = Qo[burn_in:end]
+Qo = hist_streamflow[:, "410730"]
 res_plots = []
 
 for ((res, opt), node, n_name) in zip(result, node_list, node_names)
     update_params!(node, best_candidate(res)...)
     reset!(node)
-    run_node!(node, climate)
+    @info typeof(node)
+    @time run_node!(node, climate)
 
     node_burn = node.outflow[burn_in:end]
 
@@ -68,4 +70,4 @@ combined_plot = plot(
 
 display(combined_plot)
 
-savefig("multi_model_comparison_updated.png")
+# savefig("multi_model_comparison_updated.png")
diff --git a/src/Nodes/GR4J/GR4JNode.jl b/src/Nodes/GR4J/GR4JNode.jl
index 1d4aee92..fdd8c957 100644
--- a/src/Nodes/GR4J/GR4JNode.jl
+++ b/src/Nodes/GR4J/GR4JNode.jl
@@ -46,6 +46,8 @@ abstract type GRNJNode <: NetworkNode end
 """
 GR4J Node
 
+GR4J: modèle du Génie Rural à 4 paramètres Journalier.
+
 A four-parameter model with two stores.
 
 # Parameters
@@ -287,7 +289,7 @@ end
         p_store=0.0, r_store=0.0
     )::Tuple where {F<:Float64}
 
-Generated simulated streamflow for given rainfall and potential evaporation.
+Generated simulated streamflow with GR4J for given rainfall and potential evaporation.
 
 # Parameters
 - `P` : Catchment average rainfall
diff --git a/src/Nodes/IHACRES/IHACRESExpuhNode.jl b/src/Nodes/IHACRES/_IHACRESExpuhNode.jl
similarity index 86%
rename from src/Nodes/IHACRES/IHACRESExpuhNode.jl
rename to src/Nodes/IHACRES/_IHACRESExpuhNode.jl
index 5dd84832..416194c3 100644
--- a/src/Nodes/IHACRES/IHACRESExpuhNode.jl
+++ b/src/Nodes/IHACRES/_IHACRESExpuhNode.jl
@@ -1,8 +1,14 @@
+"""
+Current disabled. 
+
+Expuh node implementation needs to be updated to use pure julia methods.
+"""
+
 using Parameters
 using ModelParameters
 
 
-Base.@kwdef mutable struct ExpuhNode{P, A<:AbstractFloat} <: IHACRESNode
+Base.@kwdef mutable struct ExpuhNode{P,A<:AbstractFloat} <: IHACRESNode
     name::String
     area::A
 
@@ -95,9 +101,9 @@ end
 
 
 function ExpuhNode(name::String, area::Float64, d::Float64, d2::Float64, e::Float64, f::Float64,
-                    tau_q::Float64, tau_s::Float64, v_s::Float64, s_coef::Float64,
-                    store::Float64, quick::Float64, slow::Float64)
-    return ExpuhNode{Param, Float64}(
+    tau_q::Float64, tau_s::Float64, v_s::Float64, s_coef::Float64,
+    store::Float64, quick::Float64, slow::Float64)
+    return ExpuhNode{Param,Float64}(
         name=name,
         area=area,
         d=d,
@@ -165,23 +171,21 @@ function run_node!(
 
     flow_res = [0.0, 0.0, 0.0]
     @ccall IHACRES.calc_flows(flow_res::Ptr{Cdouble}, prev_q::Cdouble, prev_s::Cdouble, s_node.v_s::Cdouble, e_rainfall::Cdouble,
-                              s_node.area::Cdouble, s_node.tau_q::Cdouble, s_node.tau_s::Cdouble)::Cvoid
+        s_node.area::Cdouble, s_node.tau_q::Cdouble, s_node.tau_s::Cdouble)::Cvoid
     (quick_store, slow_store, outflow) = flow_res
 
     gw_store = s_node.gw_store[end]
     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
+        routing_res::Ptr{Cdouble},
+        gw_store::Cdouble,
+        s_node.storage_coef::Cdouble,
+        inflow::Cdouble,
+        outflow::Cdouble,
+        ext::Cdouble,
+        gw_exchange::Cdouble)::Cvoid
     (gw_store, outflow) = routing_res
 
-    # level::Float64 = @ccall IHACRES.calc_ft_level(outflow::Cdouble, s_node.level_params::Ptr{Cdouble})::Cdouble
-
     update_state!(s_node, cmd, e_rainfall, et, quick_store, slow_store, outflow, gw_store)
 
     return outflow
@@ -211,19 +215,19 @@ function run_node_with_temp!(s_node::ExpuhNode, rain::Float64, temp::Float64, in
 
     flow_res = [0.0, 0.0, 0.0]
     @ccall IHACRES.calc_flows(flow_res::Ptr{Cdouble}, prev_q::Cdouble, prev_s::Cdouble, s_node.v_s::Cdouble, e_rainfall::Cdouble,
-                              s_node.area::Cdouble, s_node.tau_q::Cdouble, s_node.tau_s::Cdouble)::Cvoid
+        s_node.area::Cdouble, s_node.tau_q::Cdouble, s_node.tau_s::Cdouble)::Cvoid
     (quick_store, slow_store, outflow) = flow_res
 
     gw_store = s_node.gw_store[end]
     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
+        routing_res::Ptr{Cdouble},
+        gw_store::Cdouble,
+        s_node.storage_coef::Cdouble,
+        inflow::Cdouble,
+        outflow::Cdouble,
+        ext::Cdouble,
+        gw_exchange::Cdouble)::Cvoid
     (gw_store, outflow) = routing_res
 
     # level::Float64 = @ccall IHACRES.calc_ft_level(outflow::Cdouble, s_node.level_params::Ptr{Cdouble})::Cdouble
@@ -235,7 +239,7 @@ end
 
 
 function update_params!(node::ExpuhNode, d::Float64, d2::Float64, e::Float64, f::Float64,
-                        tau_q::Float64, tau_s::Float64, v_s::Float64, s_coef::Float64)::Nothing
+    tau_q::Float64, tau_s::Float64, v_s::Float64, s_coef::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)
diff --git a/src/Streamfall.jl b/src/Streamfall.jl
index e1a562fc..f4aea096 100644
--- a/src/Streamfall.jl
+++ b/src/Streamfall.jl
@@ -6,29 +6,13 @@ using Graphs, MetaGraphs, Distributed, DataFrames
 
 const MODPATH = @__DIR__
 
-if Sys.iswindows()
-    libext = ".dll"
-elseif Sys.islinux()
-    libext = ".so"
-elseif Sys.isapple()
-    libext = ".dynlib"
-else
-    throw(DomainError("Unsupported platform"))
-end
-
-# Can't use string, DLL location has to be a const
-# (which makes sense but still, many hours wasted!)
-# https://github.com/JuliaLang/julia/issues/29602
-const IHACRES = joinpath(MODPATH, "../deps", "usr", "lib", "ihacres$(libext)")
-
-
 include("Network.jl")
 include("Nodes/Node.jl")
 include("Climate.jl")
 include("metrics.jl")
 include("calibration.jl")
 include("Nodes/IHACRES/IHACRESNode.jl")
-include("Nodes/IHACRES/IHACRESExpuhNode.jl")
+# include("Nodes/IHACRES/IHACRESExpuhNode.jl")
 include("Nodes/GR4J/GR4JNode.jl")
 include("Nodes/HyMod/HyModNode.jl")
 include("Nodes/SYMHYD/SYMHYDNode.jl")
@@ -37,9 +21,9 @@ include("Nodes/Ensembles/EnsembleNode.jl")
 
 
 function timestep_value(ts::Int, gauge_id::String, col_partial::String, dataset::DataFrame)::Float64
-    target_col = filter(x -> occursin(gauge_id, x)
-                                & occursin(col_partial, x),
-                                names(dataset)
+    target_col = filter(
+        x -> occursin(gauge_id, x) & occursin(col_partial, x),
+        names(dataset)
     )
 
     amount = 0.0
@@ -107,7 +91,7 @@ julia> climate, streamflow = align_time_frame(climate, streamflow)
 """
 function align_time_frame(timeseries::T...) where {T<:DataFrame}
     min_date, max_date = find_common_timeframe(timeseries...)
-    modded = [t[min_date .<= t.Date .<= max_date, :] for t in timeseries]
+    modded = [t[min_date.<=t.Date.<=max_date, :] for t in timeseries]
 
     for t in modded
         ts_diff = diff(t.Date)
@@ -137,7 +121,7 @@ function run_basin!(
     _, outlets = find_inlets_and_outlets(sn)
     @inbounds for outlet_id in outlets
         run_node!(sn, outlet_id, climate;
-                  inflow=inflow, extraction=extraction, exchange=exchange)
+            inflow=inflow, extraction=extraction, exchange=exchange)
     end
 end
 
@@ -294,9 +278,9 @@ include("plotting.jl")
 
 # Nodes
 export NetworkNode, GenericNode, GenericDirectNode
-export IHACRES, IHACRESNode, IHACRESBilinearNode, ExpuhNode, DamNode, Climate
+export IHACRES, IHACRESNode, IHACRESBilinearNode, DamNode, Climate
 export create_node, GR4JNode, HyModNode, SimpleHyModNode, SYMHYDNode
-export EnsembleNode, WeightedEnsembleNode
+export EnsembleNode, WeightedEnsembleNode, GREnsembleNode
 
 # Network
 export find_inlets_and_outlets, inlets, outlets
diff --git a/src/metrics.jl b/src/metrics.jl
index dad14ec0..f46fb697 100644
--- a/src/metrics.jl
+++ b/src/metrics.jl
@@ -139,7 +139,10 @@ end
 """
     NSE(obs, sim)
 
-The Nash-Sutcliffe Efficiency score
+The Nash-Sutcliffe Efficiency score.
+
+Ranges from 1.0 to -∞, where values below zero indicate the mean of
+observations would outperform a model.
 """
 function NSE(obs, sim)
     return 1.0 - sum((obs .- sim) .^ 2) / sum((obs .- mean(obs)) .^ 2)
@@ -164,7 +167,7 @@ end
 """
     RMSE(obs, sim)
 
-Root Mean Square Error
+Root Mean Square Error.
 """
 function RMSE(obs, sim)
     return (sum((sim .- obs) .^ 2) / length(sim))^0.5
@@ -186,7 +189,7 @@ end
 """
     ADJ_R2(obs, sim, p::Int64)::Float64
 
-Adjusted R²
+Adjusted R².
 
 # Arguments
 - `obs::Vector` : observations
@@ -204,7 +207,7 @@ end
 """
     MAE(obs, sim)
 
-Mean Absolute Error
+Mean Absolute Error.
 """
 MAE(obs, sim) = mean(abs.(sim .- obs))
 
@@ -212,7 +215,7 @@ MAE(obs, sim) = mean(abs.(sim .- obs))
 """
     ME(obs, sim)
 
-Mean Error
+Mean Error.
 """
 ME(obs, sim) = mean(sim .- obs)
 
@@ -471,7 +474,7 @@ function mKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
     # Timing (Pearson's correlation)
     r = Statistics.cor(obs, sim)
     if isnan(r)
-        # can happen if obs or sim is a constant (std of 0)
+        # Can happen if obs or sim is a constant (std of 0)
         r = all(obs .== sim) ? 1.0 : 0.0
     end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index ff056b4d..5695ddd8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,7 +10,7 @@ TEST_DIR = @__DIR__
     ihacres = create_node(IHACRESBilinearNode, "IHACRES", 100.0)
 
     # Expuh form does not yet support time stepping
-    expuh = create_node(ExpuhNode, "Expuh", 100.0)
+    # expuh = create_node(ExpuhNode, "Expuh", 100.0)
 
     gr4j = create_node(GR4JNode, "GR4J", 100.0)
     hymod = create_node(SimpleHyModNode, "HyMod", 100.0)
@@ -22,8 +22,8 @@ TEST_DIR = @__DIR__
         Streamfall.prep_state!(ihacres, 1)
         @test Streamfall.run_timestep!(ihacres, 6.0, 3.0, 1) isa AbstractFloat
 
-        Streamfall.prep_state!(expuh, 1)
-        @test Streamfall.run_node!(expuh, 6.0, 3.0, 50.0, 10.0, 5.0) isa AbstractFloat
+        # Streamfall.prep_state!(expuh, 1)
+        # @test Streamfall.run_node!(expuh, 6.0, 3.0, 50.0, 10.0, 5.0) isa AbstractFloat
 
         Streamfall.prep_state!(gr4j, 1)
         @test Streamfall.run_timestep!(gr4j, 6.0, 3.0, 1) isa AbstractFloat
@@ -86,15 +86,13 @@ end
     params = (214.6561105573191, 76.6251447, 200.0, 2.0, 0.727)
     current_store, rain, d, d2, alpha = params
 
-    interim_results = [0.0, 0.0, 0.0]
-    @ccall IHACRES.calc_ft_interim_cmd(interim_results::Ptr{Cdouble},
-                                       current_store::Cdouble,
-                                       rain::Cdouble,
-                                       d::Cdouble,
-                                       d2::Cdouble,
-                                       alpha::Cdouble)::Cvoid
-
-    (mf, e_rainfall, recharge) = interim_results
+    (mf, e_rainfall, recharge) = Streamfall.calc_ft_interim_cmd(
+        current_store,
+        rain,
+        d,
+        d2,
+        alpha
+    )
 
     @test !isnan(mf)
     @test !isnan(e_rainfall)
@@ -108,73 +106,73 @@ end
     recharge = 3.84930005080411E-06
     rain = 0.0000188
 
-    n_cmd = @ccall IHACRES.calc_cmd(cmd::Cdouble, rain::Cdouble, et::Cdouble, e_rain::Cdouble, recharge::Cdouble)::Float64
+    n_cmd = Streamfall.calc_cmd(cmd, rain, et, e_rain, recharge)
 
     @test isapprox(n_cmd, 106.22, atol=0.001)
 end
 
 
-@testset "IHACRES calculations" begin
-    area = 1985.73
-    a = 54.352
-    b = 0.187
-    e_rain = 3.421537294474909e-6
-    recharge = 3.2121031313153022e-6
-
-    prev_quick = 100.0
-    prev_slow = 100.0
-
-    flow_results = [0.0, 0.0, 0.0]
-    @ccall IHACRES.calc_ft_flows(
-        flow_results::Ptr{Cdouble},
-        prev_quick::Cdouble,
-        prev_slow::Cdouble,
-        e_rain::Cdouble,
-        recharge::Cdouble,
-        area::Cdouble,
-        a::Cdouble,
-        b::Cdouble
-    )::Cvoid
-
-    quickflow = (prev_quick + (e_rain * area))
-    α = exp(-a)
-    quickflow = α * quickflow
-
-    @test flow_results[1] == quickflow
-
-    slow_store = prev_slow + (recharge * area)
-    α = exp(-b)
-    slow_store = α * slow_store
-    @test flow_results[2] == slow_store
-
-    e_rain = 0.0
-    recharge = 0.0
-
-    prev_quick = 3.3317177943791187
-    prev_slow = 144.32012122323678
-
-    flow_results = [0.0, 0.0, 0.0]
-    @ccall IHACRES.calc_ft_flows(
-        flow_results::Ptr{Cdouble},
-        prev_quick::Cdouble,
-        prev_slow::Cdouble,
-        e_rain::Cdouble,
-        recharge::Cdouble,
-        area::Cdouble,
-        a::Cdouble,
-        b::Cdouble
-    )::Cvoid
-
-    quickflow = (prev_quick + (e_rain * area))
-    α = exp(-a)
-    @test flow_results[1] == (α * quickflow)
-
-    slow_store = prev_slow + (recharge * area)
-    α = exp(-b)
-    β = (1.0 - α) * slow_store
-    slow_store = α * slow_store
-    @test flow_results[2] == slow_store
-end
+# @testset "IHACRES calculations" begin
+#     area = 1985.73
+#     a = 54.352
+#     b = 0.187
+#     e_rain = 3.421537294474909e-6
+#     recharge = 3.2121031313153022e-6
+
+#     prev_quick = 100.0
+#     prev_slow = 100.0
+
+#     flow_results = [0.0, 0.0, 0.0]
+#     @ccall IHACRES.calc_ft_flows(
+#         flow_results::Ptr{Cdouble},
+#         prev_quick::Cdouble,
+#         prev_slow::Cdouble,
+#         e_rain::Cdouble,
+#         recharge::Cdouble,
+#         area::Cdouble,
+#         a::Cdouble,
+#         b::Cdouble
+#     )::Cvoid
+
+#     quickflow = (prev_quick + (e_rain * area))
+#     α = exp(-a)
+#     quickflow = α * quickflow
+
+#     @test flow_results[1] == quickflow
+
+#     slow_store = prev_slow + (recharge * area)
+#     α = exp(-b)
+#     slow_store = α * slow_store
+#     @test flow_results[2] == slow_store
+
+#     e_rain = 0.0
+#     recharge = 0.0
+
+#     prev_quick = 3.3317177943791187
+#     prev_slow = 144.32012122323678
+
+#     flow_results = [0.0, 0.0, 0.0]
+#     @ccall IHACRES.calc_ft_flows(
+#         flow_results::Ptr{Cdouble},
+#         prev_quick::Cdouble,
+#         prev_slow::Cdouble,
+#         e_rain::Cdouble,
+#         recharge::Cdouble,
+#         area::Cdouble,
+#         a::Cdouble,
+#         b::Cdouble
+#     )::Cvoid
+
+#     quickflow = (prev_quick + (e_rain * area))
+#     α = exp(-a)
+#     @test flow_results[1] == (α * quickflow)
+
+#     slow_store = prev_slow + (recharge * area)
+#     α = exp(-b)
+#     β = (1.0 - α) * slow_store
+#     slow_store = α * slow_store
+#     @test flow_results[2] == slow_store
+# end
 
 include("test_metrics.jl")
 include("test_data_op.jl")