diff --git a/README.md b/README.md
index 81c38113..82819cde 100644
--- a/README.md
+++ b/README.md
@@ -62,13 +62,6 @@ julia>] test
 
 ## Usage
 
-Streamfall is currently unregistered but can be added to a Julia environment directly from
-the Package manager:
-
-```bash
-julia>] add https://github.com/ConnectedSystems/Streamfall.jl#main
-```
-
 The examples below use data from the CAMEL-AUS dataset, available here:
 
 > Fowler, K. J. A., Acharya, S. C., Addor, N., Chou, C., and Peel, M. C.: CAMELS-AUS: hydrometeorological time   series and landscape attributes for 222 catchments in Australia, Earth Syst. Sci. Data, 13, 3847–3867, https://doi.org/10.5194/essd-13-3847-2021, 2021.
diff --git a/src/Climate.jl b/src/Climate.jl
index b2c50c37..c37aa973 100644
--- a/src/Climate.jl
+++ b/src/Climate.jl
@@ -19,6 +19,10 @@ function Climate(file_path::String, p_id, et_id; t_id="_T")
 end
 
 """
+    extract_flow(
+        data::DataFrame, gauge_id::String; suffix::String="_Q"
+    )::DataFrame
+
 Extract streamflow data from file.
 
 Streamflow (Q) column is identified the Gauge ID.
@@ -35,13 +39,13 @@ e.g., ("000001_Q")
 DataFrame of observations for selected gauge.
 """
 @inline function extract_flow(
-    data::DataFrame, gauge_id::String, suffix::String="_Q"
+    data::DataFrame, gauge_id::String; suffix::String="_Q"
 )::DataFrame
     target = data[:, ["Date", gauge_id * suffix]]
     try
-        target[!, gauge_id * suffix] = convert.(Float64, target[!, gauge_id * suffix])
+        target[!, gauge_id*suffix] = convert.(Float64, target[!, gauge_id*suffix])
     catch
-        target[!, gauge_id * suffix] = convert.(Union{Float64,Missing}, target[!, gauge_id * suffix])
+        target[!, gauge_id*suffix] = convert.(Union{Float64,Missing}, target[!, gauge_id*suffix])
     end
 
     rename!(target, gauge_id * suffix => gauge_id)
@@ -78,8 +82,9 @@ Extract rainfall data for a given node.
 function rainfall_data(node::NetworkNode, climate::Climate)::DataFrame
     data = climate.climate_data
     rain_col = filter(x -> occursin(node.name, x)
-                            & occursin(climate.rainfall_id, x),
-                            names(data))[1]
+                           &
+                           occursin(climate.rainfall_id, x),
+        names(data))[1]
 
     return data[:, rain_col]
 end
@@ -107,15 +112,18 @@ function climate_values(node::NetworkNode, climate::Climate, timestep::Int)
 
     # TODO : Catch instances where data is not found (raises BoundsError)
     rain_col = filter(x -> occursin(node_name, x)
-                        & occursin(climate.rainfall_id, x),
-                        names(data))[1]
+                           &
+                           occursin(climate.rainfall_id, x),
+        names(data))[1]
     et_col = filter(x -> occursin(node_name, x)
-                        & occursin(climate.et_id, x),
-                        names(data))[1]
+                         &
+                         occursin(climate.et_id, x),
+        names(data))[1]
     t_col = try
         filter(x -> occursin(node_name, x)
-                            & occursin(climate.t_id, x),
-                            names(data))[1]
+                    &
+                    occursin(climate.t_id, x),
+            names(data))[1]
     catch err
         if !(err isa BoundsError)
             rethrow(err)
@@ -148,11 +156,13 @@ function climate_values(node::NetworkNode, climate::Climate)
 
     # TODO : Catch instances where data is not found (raises BoundsError)
     rain_col = filter(x -> occursin(node.name, x)
-                            & occursin(climate.rainfall_id, x),
-                            names(data))[1]
+                           &
+                           occursin(climate.rainfall_id, x),
+        names(data))[1]
     et_col = filter(x -> occursin(node.name, x)
-                            & occursin(climate.et_id, x),
-                            names(data))[1]
+                         &
+                         occursin(climate.et_id, x),
+        names(data))[1]
 
     if isempty(rain_col) | isempty(et_col)
         throw(ArgumentError("No climate data found for $(node.name) at time step: $(timestep)"))
diff --git a/src/Network.jl b/src/Network.jl
index 3fd0931f..2d47f982 100644
--- a/src/Network.jl
+++ b/src/Network.jl
@@ -41,7 +41,11 @@ function get_prop(sn::StreamfallNetwork, nid::Int64, prop::Symbol)::Any
 end
 
 
-"""Determine a node's connection"""
+"""
+    in_or_out(g, v)
+
+Determine a node's connection
+"""
 function in_or_out(g, v)
     ins = length(inneighbors(g, v))
     outs = length(outneighbors(g, v))
@@ -58,7 +62,11 @@ function in_or_out(g, v)
 end
 
 
-"""Find all inlets and outlets in a network."""
+"""
+    find_inlets_and_outlets(sn::StreamfallNetwork)::Tuple
+
+Find all inlets and outlets in a network.
+"""
 function find_inlets_and_outlets(sn::StreamfallNetwork)::Tuple
     g = sn.mg
     vs = vertices(g)
@@ -161,9 +169,14 @@ function create_node(mg::MetaDiGraph, node_name::String, details::AbstractDict,
             func = run_node!
         end
 
-        set_props!(mg, nid, Dict(:name=>node_name,
-                                 :node=>n,
-                                 :nfunc=>func))
+        set_props!(
+            mg, nid,
+            Dict(
+                :name => node_name,
+                :node => n,
+                :nfunc => func
+            )
+        )
 
         this_id = nid
     else
@@ -307,7 +320,7 @@ function Base.show(io::IO, sn::StreamfallNetwork)
 
     show_verts = vs
     if length(vs) > 4
-        show_verts = [1, 2, nothing, length(vs)-1, length(vs)]
+        show_verts = [1, 2, nothing, length(vs) - 1, length(vs)]
     end
 
     for nid in show_verts
@@ -328,7 +341,7 @@ end
 Simple plot of stream network.
 """
 function plot_network(sn::StreamfallNetwork; as_html=false)
-    node_labels = ["$(sn[i].name)\n"*string(nameof(typeof(sn[i]))) for i in vertices(sn.mg)]
+    node_labels = ["$(sn[i].name)\n" * string(nameof(typeof(sn[i]))) for i in vertices(sn.mg)]
 
     if as_html
         plot_func = gplothtml
@@ -355,7 +368,7 @@ function Base.iterate(sn::StreamfallNetwork, state=1)
         return nothing
     end
 
-    return (sn[state], state+1)
+    return (sn[state], state + 1)
 end
 
 
diff --git a/src/Nodes/DamNode.jl b/src/Nodes/DamNode.jl
index d2021ae6..0518316b 100644
--- a/src/Nodes/DamNode.jl
+++ b/src/Nodes/DamNode.jl
@@ -26,7 +26,7 @@ function c_dam_outflow(discharge, irrigation_extraction)
 end
 
 
-Base.@kwdef mutable struct DamNode{P, A<:AbstractFloat} <: NetworkNode
+Base.@kwdef mutable struct DamNode{P,A<:AbstractFloat} <: NetworkNode
     name::String
     area::A
 
@@ -74,12 +74,12 @@ function DamNode(
     calc_dam_outflow::Function
 ) where {F<:Float64}
     return DamNode(name, area, max_storage, storage_coef,
-                   calc_dam_level, calc_dam_area, calc_dam_discharge, calc_dam_outflow,
-                   F[initial_storage], F[], F[], F[], F[], F[], F[], F[])
+        calc_dam_level, calc_dam_area, calc_dam_discharge, calc_dam_outflow,
+        F[initial_storage], F[], F[], F[], F[], F[], F[], F[])
 end
 
 """
-    DamNode(name::String, spec::Dict)
+    DamNode(name::String, spec::AbstractDict)
 
 Create DamNode from a given specification.
 """
@@ -136,7 +136,7 @@ function storage(node::DamNode)
 end
 
 function prep_state!(node::DamNode, timesteps::Int64)
-    resize!(node.storage, timesteps+1)
+    resize!(node.storage, timesteps + 1)
     node.storage[2:end] .= 0.0
 
     node.effective_rainfall = zeros(timesteps)
@@ -204,13 +204,13 @@ end
 
 
 function run_node!(node::DamNode, climate::Climate;
-                   inflow=nothing, extraction=nothing, exchange=nothing)
+    inflow=nothing, extraction=nothing, exchange=nothing)
     timesteps = sim_length(climate)
     prep_state!(node, timesteps)
 
     for ts in 1:timesteps
         run_node!(node, climate, ts;
-                  inflow=inflow, extraction=extraction, exchange=exchange)
+            inflow=inflow, extraction=extraction, exchange=exchange)
     end
 
     return nothing
@@ -289,7 +289,7 @@ function run_node!(
     discharge = node.calc_dam_discharge(volume, node.max_storage)
 
     updated_store = update_volume(volume, inflow, gw_flux, rain, et,
-                                  dam_area, extractions, discharge, node.max_storage)
+        dam_area, extractions, discharge, node.max_storage)
     outflow = node.calc_dam_outflow(discharge, extractions)
 
     update_state!(node, ts, updated_store, rain, et, dam_area, discharge, outflow)
diff --git a/src/Nodes/GR4J/GR4JNode.jl b/src/Nodes/GR4J/GR4JNode.jl
index 2df0d5d0..1d4aee92 100644
--- a/src/Nodes/GR4J/GR4JNode.jl
+++ b/src/Nodes/GR4J/GR4JNode.jl
@@ -49,10 +49,10 @@ GR4J Node
 A four-parameter model with two stores.
 
 # Parameters
-- x1 : maximum capacity of the production store (mm) (> 0)
-- x2 : groundwater exchange coefficient (mm) (value < and > 0 possible)
-- x3 : one day ahead maximum capacity of the routing store (mm, > 0)
-- x4 : time base of unit hydrograph UH1 (days, > 0.5)
+- `x1` : maximum capacity of the production store (mm) (> 0)
+- `x2` : groundwater exchange coefficient (mm) (value < and > 0 possible)
+- `x3` : one day ahead maximum capacity of the routing store (mm, > 0)
+- `x4` : time base of unit hydrograph UH1 (days, > 0.5)
 
 # References
 1. Perrin, C., Michel, C., Andréassian, V., 2003.
@@ -64,20 +64,19 @@ A four-parameter model with two stores.
     Python GR4J
     https://github.com/amacd31/gr4j
 """
-Base.@kwdef mutable struct GR4JNode{P, A<:AbstractFloat} <: GRNJNode
+Base.@kwdef mutable struct GR4JNode{P,A<:AbstractFloat} <: GRNJNode
     const name::String
     const area::A
 
-    # parameters
-    X1::P = Param(350.0, bounds=(1.0, 1500.0))
-    X2::P = Param(0.0, bounds=(-10.0, 5.0))
-    X3::P = Param(40.0, bounds=(1.0, 500.0))
-    X4::P = Param(0.5, bounds=(0.5, 10.0))
-
+    # Parameters
     # x1 : maximum capacity of the production store (mm) (> 0)
     # x2 : groundwater exchange coefficient (mm) (value < and > 0 possible)
     # x3 : one day ahead maximum capacity of the routing store (mm, > 0)
     # x4 : time base of unit hydrograph UH1 (days, > 0.5)
+    X1::P = Param(350.0, bounds=(1.0, 1500.0))
+    X2::P = Param(0.0, bounds=(-10.0, 5.0))
+    X3::P = Param(40.0, bounds=(1.0, 500.0))
+    X4::P = Param(0.5, bounds=(0.5, 10.0))
 
     # stores
     p_store::Vector{A} = [0.0]
@@ -96,15 +95,15 @@ end
 
 
 function prep_state!(node::GR4JNode, timesteps::Int64)
-    resize!(node.p_store, timesteps+1)
-    resize!(node.r_store, timesteps+1)
+    resize!(node.p_store, timesteps + 1)
+    resize!(node.r_store, timesteps + 1)
     node.p_store[2:end] .= 0.0
     node.r_store[2:end] .= 0.0
 
     # Prep cache
     X4 = node.X4.val
     nUH1 = Int(ceil(X4))
-    nUH2 = Int(ceil(2.0*X4))
+    nUH2 = Int(ceil(2.0 * X4))
     cUH1, cUH2 = fill(0.0, nUH1), fill(0.0, nUH2)
     node.UH1 = fill(cUH1, timesteps)
     node.UH2 = fill(cUH2, timesteps)
@@ -114,7 +113,6 @@ function prep_state!(node::GR4JNode, timesteps::Int64)
     node.outflow = fill(0.0, timesteps)
 end
 
-
 function GR4JNode(name::String, spec::AbstractDict)
     n = create_node(GR4JNode, name, spec["area"])
     node_params = spec["parameters"]
@@ -173,7 +171,11 @@ function run_timestep!(
 )
     uh1_cache = node.uh1_ordinates
     uh2_cache = node.uh2_ordinates
-    res = run_gr4j(rain, et, node.X1.val, node.X2.val, node.X3.val, node.X4.val, node.area, node.UH1[ts], node.UH2[ts], uh1_cache, uh2_cache, node.p_store[ts], node.r_store[ts])
+    res = run_gr4j(
+        rain, et, node.X1.val, node.X2.val, node.X3.val, node.X4.val, node.area,
+        node.UH1[ts], node.UH2[ts], uh1_cache, uh2_cache;
+        p_store=node.p_store[ts], r_store=node.r_store[ts]
+    )
     Q, p_s, r_s, UH1, UH2 = res
 
     node_name = node.name
@@ -190,13 +192,20 @@ function run_timestep!(
     return Q
 end
 
+"""
+    update_state!(node::GR4JNode, ps, rs, q, UH1, UH2)::Nothing
+    update_state!(node::GR4JNode, ts::Int64, ps, rs, q, UH1, UH2)::Nothing
 
-function update_state!(node::GR4JNode, ps, rs, q, UH1, UH2)
+Update GR4J node state.
+"""
+function update_state!(node::GR4JNode, ps, rs, q, UH1, UH2)::Nothing
     append!(node.p_store, ps)
     append!(node.r_store, rs)
     append!(node.outflow, q)
     push!(node.UH1, UH1)
     push!(node.UH2, UH2)
+
+    return nothing
 end
 function update_state!(node::GR4JNode, ts::Int64, ps, rs, q, UH1, UH2)::Nothing
     node.p_store[ts+1] = ps
@@ -230,7 +239,7 @@ end
 Reset node. Clears all states back to their initial values.
 """
 function reset!(node::GR4JNode)::Nothing
-    # stores
+    # Stores
     node.p_store = [node.p_store[1]]
     node.r_store = [node.r_store[1]]
 
@@ -248,17 +257,17 @@ end
 
 Determine unit hydrograph ordinates.
 """
-function s_curve(t::F, x4::F; uh2::Bool = false)::F where {F<:Float64}
+function s_curve(t::F, x4::F; uh2::Bool=false)::F where {F<:Float64}
     if t <= 0.0
         return 0.0
     end
 
     ordinate::F = 0.0
     if t < x4
-        ordinate = (t/x4)^2.5
+        ordinate = (t / x4)^2.5
     else
-        if uh2 && (t < 2*x4)
-            ordinate = 1.0 - 0.5*(2 - t/x4)^2.5
+        if uh2 && (t < 2 * x4)
+            ordinate = 1.0 - 0.5 * (2 - t / x4)^2.5
         else
             # t >= x4 if uh1, or
             # t >= 2*x4 if uh2
@@ -271,56 +280,77 @@ end
 
 
 """
-    run_gr4j(P::Float64, E::Float64,
-             X1::Float64, X2::Float64, X3::Float64, X4::Float64, area::Float64,
-             p_store::Float64=0.0, r_store::Float64=0.0)::Tuple
+    run_gr4j(
+        P::F, E::F, X1::F, X2::F, X3::F, X4::F, area::F,
+        UH1::Vector{Float64}, UH2::Vector{Float64},
+        uh1_ordinates::Vector{Float64}, uh2_ordinates::Vector{Float64};
+        p_store=0.0, r_store=0.0
+    )::Tuple where {F<:Float64}
 
 Generated simulated streamflow for given rainfall and potential evaporation.
 
 # Parameters
-- P : Catchment average rainfall
-- E : Catchment average potential evapotranspiration
-- X1 - X4 : X parameters
-- area : Catchment area
-- p_store : Initial production store
-- r_store : Initial state store
+- `P` : Catchment average rainfall
+- `E` : Catchment average potential evapotranspiration
+- `X1` : Maximum capacity of production store (in mm; > 0)
+- `X2` : Groundwater exchange coefficient (in mm; value < and > 0 possible)
+- `X3` : Maximum capacity of routing store (in mm; > 0)
+- `X4` : Time base of the unit hydrograph (in days, > 0.5)
+- `area` : Catchment area
+- `UH1` : Quickflow store
+- `UH2` : Baseflow store
+- `uh1_ordinates` : The proportion of rainfall converted to quickflow for each timestep
+- `uh2_ordinates` : The proportion of rainfall converted to slowflow for each timestep
+- `p_store` : Initial production store
+- `r_store` : Initial state store
 
 # Returns
-- tuple of simulated outflow [ML/day], and intermediate states: p_store, r_store, UH1, UH2
+Tuple of:
+- Simulated outflow [ML/day]
+- intermediate states:
+  - p_store (initial production / percolation)
+  - r_store (initial state)
+  - UH1 (Quickflow)
+  - UH2 (Slowflow)
 """
-function run_gr4j(P::F, E::F, X1::F, X2::F, X3::F, X4::F, area::F, UH1::Vector{Float64}, UH2::Vector{Float64}, uh1_ordinates::Vector{Float64}, uh2_ordinates::Vector{Float64}, p_store=0.0, r_store=0.0)::Tuple where {F<:Float64}
+function run_gr4j(
+    P::F, E::F, X1::F, X2::F, X3::F, X4::F, area::F,
+    UH1::Vector{F}, UH2::Vector{F},
+    uh1_ordinates::Vector{F}, uh2_ordinates::Vector{F};
+    p_store=0.0, r_store=0.0
+)::Tuple where {F<:Float64}
     nUH1::Int64 = Int(ceil(X4))
-    nUH2::Int64 = Int(ceil(2.0*X4))
+    nUH2::Int64 = Int(ceil(2.0 * X4))
 
     @inbounds for t in 2:(nUH1+1)
         t_f = Float64(t)
-        uh1_ordinates[t - 1] = s_curve(t_f, X4) - s_curve(t_f-1.0, X4)
+        uh1_ordinates[t-1] = s_curve(t_f, X4) - s_curve(t_f - 1.0, X4)
     end
 
     @inbounds for t in 2:(nUH2+1)
         t_f = Float64(t)
-        uh2_ordinates[t - 1] = s_curve(t_f, X4, uh2=true) - s_curve(t_f-1.0, X4, uh2=true)
+        uh2_ordinates[t-1] = s_curve(t_f, X4, uh2=true) - s_curve(t_f - 1.0, X4, uh2=true)
     end
 
     Q::F = 0.0
     if P > E
         net_evap = 0.0
-        scaled_net_precip = min((P - E)/X1, 13.0)
+        scaled_net_precip = min((P - E) / X1, 13.0)
         tanh_scaled_net_precip = tanh(scaled_net_precip)
 
-        tmp_a = (p_store/X1)^2
+        tmp_a = (p_store / X1)^2
         numer = (X1 * (1.0 - tmp_a) * tanh_scaled_net_precip)
-        denom = (1.0 + p_store/X1 * tanh_scaled_net_precip)
+        denom = (1.0 + p_store / X1 * tanh_scaled_net_precip)
         reservoir_production = numer / denom
 
         routed_volume = P - E - reservoir_production
     else
-        scaled_net_evap = min((E - P)/X1, 13.0)
+        scaled_net_evap = min((E - P) / X1, 13.0)
         tanh_scaled_net_evap = tanh(scaled_net_evap)
 
-        ps_div_x1 = (2.0 - p_store/X1) * tanh_scaled_net_evap
+        ps_div_x1 = (2.0 - p_store / X1) * tanh_scaled_net_evap
         net_evap = (p_store * (ps_div_x1) /
-                    (1.0 + (1.0 - p_store/X1) * tanh_scaled_net_evap))
+                    (1.0 + (1.0 - p_store / X1) * tanh_scaled_net_evap))
 
         reservoir_production = 0.0
         routed_volume = 0.0
@@ -332,15 +362,15 @@ function run_gr4j(P::F, E::F, X1::F, X2::F, X3::F, X4::F, area::F, UH1::Vector{F
     tmp_b::F = (1 + tmp_a)^0.25
     percolation = p_store / tmp_b
 
-    routed_volume = routed_volume + (p_store-percolation)
+    routed_volume = routed_volume + (p_store - percolation)
 
     @inbounds for i in 1:nUH1-1
-        UH1[i] = UH1[i+1] + uh1_ordinates[i]*routed_volume
+        UH1[i] = UH1[i+1] + uh1_ordinates[i] * routed_volume
     end
     UH1[end] = uh1_ordinates[end] * routed_volume
 
     @inbounds for j in 1:nUH2-1
-        UH2[j] = UH2[j+1] + uh2_ordinates[j]*routed_volume
+        UH2[j] = UH2[j+1] + uh2_ordinates[j] * routed_volume
     end
     UH2[end] = uh2_ordinates[end] * routed_volume
 
@@ -353,7 +383,7 @@ function run_gr4j(P::F, E::F, X1::F, X2::F, X3::F, X4::F, area::F, UH1::Vector{F
     R2 = r_store / tmp_b
     QR = r_store - R2
     r_store = R2
-    QD = max(0.0, UH2[1]*0.1+groundwater_exchange)
+    QD = max(0.0, UH2[1] * 0.1 + groundwater_exchange)
 
     Q = (QR + QD) * area
 
diff --git a/src/Nodes/HyMod/HyModNode.jl b/src/Nodes/HyMod/HyModNode.jl
index cd686644..30c0fae6 100644
--- a/src/Nodes/HyMod/HyModNode.jl
+++ b/src/Nodes/HyMod/HyModNode.jl
@@ -20,7 +20,7 @@ Adapted with kind permission from: https://github.com/jdherman/GRA-2020-SALib
     tion of hydrological models, Hydrol. Earth Syst. Sci., 5, 13-26,
     https://doi.org/10.5194/hess-5-13-2001.
 """
-Base.@kwdef mutable struct SimpleHyModNode{P, A<:AbstractFloat} <: HyModNode
+Base.@kwdef mutable struct SimpleHyModNode{P,A<:AbstractFloat} <: HyModNode
     const name::String
     const area::A
 
@@ -56,9 +56,16 @@ function SimpleHyModNode(name::String, spec::AbstractDict)
     return n
 end
 
-
-function SimpleHyModNode(name::String, area::Float64, sm_max::Float64, B::Float64,
-                         alpha::Float64, Kf::Float64, Ks::Float64)
+"""
+    SimpleHyModNode(
+        name::String, area::Float64, sm_max::Float64, B::Float64,
+        alpha::Float64, Kf::Float64, Ks::Float64
+    )
+"""
+function SimpleHyModNode(
+    name::String, area::Float64, sm_max::Float64, B::Float64,
+    alpha::Float64, Kf::Float64, Ks::Float64
+)
     n = create_node(SimpleHyModNode, name, area)
     update_params!(n, sm_max, B, alpha, Kf, Ks)
 
@@ -75,11 +82,11 @@ end
 
 
 function prep_state!(node::SimpleHyModNode, sim_length::Int64)::Nothing
-    resize!(node.Sm, sim_length+1)
-    resize!(node.Sf1, sim_length+1)
-    resize!(node.Sf2, sim_length+1)
-    resize!(node.Sf3, sim_length+1)
-    resize!(node.Ss1, sim_length+1)
+    resize!(node.Sm, sim_length + 1)
+    resize!(node.Sf1, sim_length + 1)
+    resize!(node.Sf2, sim_length + 1)
+    resize!(node.Sf3, sim_length + 1)
+    resize!(node.Ss1, sim_length + 1)
 
     node.Sm[2:end] .= 0.0
     node.Sf1[2:end] .= 0.0
@@ -94,8 +101,10 @@ function prep_state!(node::SimpleHyModNode, sim_length::Int64)::Nothing
 end
 
 """
-    run_timestep!(node::HyModNode, climate::Climate, timestep::Int,
-                  inflow::Float64, extraction::Float64, exchange::Float64)
+    run_timestep!(
+        node::SimpleHyModNode, climate::Climate, timestep::Int;
+        inflow=nothing, extraction=nothing, exchange=nothing
+    )::Float64
 
 Run given HyMod node for a time step.
 """
@@ -136,20 +145,20 @@ function run_hymod!(node::SimpleHyModNode, ts::Int64, P::F, PET::F, inflow::F, e
     Ks::F = node.Ks.val::F
 
     # Calculate fluxes
-    Peff::F = P*(1.0 - max(1.0 - Sm/Sm_max, 0.0)^B)  # PDM model Moore 1985
-    ETa::F = min(PET*(Sm/Sm_max), Sm)
+    Peff::F = P * (1.0 - max(1.0 - Sm / Sm_max, 0.0)^B)  # PDM model Moore 1985
+    ETa::F = min(PET * (Sm / Sm_max), Sm)
 
-    Qf1::F = Kf*Sf1
-    Qf2::F = Kf*Sf2
-    Qf3::F = Kf*Sf3
-    Qs1::F = Ks*Ss1
+    Qf1::F = Kf * Sf1
+    Qf2::F = Kf * Sf2
+    Qf3::F = Kf * Sf3
+    Qs1::F = Ks * Ss1
 
     # update state variables
     Sm_t1::F = max(Sm + P - Peff - ETa, 0.0)
-    Sf1_t1::F = Sf1 + alpha*Peff - Qf1
+    Sf1_t1::F = Sf1 + alpha * Peff - Qf1
     Sf2_t1::F = Sf2 + Qf1 - Qf2
     Sf3_t1::F = Sf3 + Qf2 - Qf3
-    Ss1_t1::F = Ss1 + (1-alpha)*Peff - Qs1
+    Ss1_t1::F = Ss1 + (1 - alpha) * Peff - Qs1
 
     Qtmp::F = (Qs1 + Qf3) * node.area
     Q_t::F = inflow + Qtmp - ext + flux
@@ -171,17 +180,6 @@ function update_params!(node::HyModNode, Sm_max::F, B::F, alpha::F, Kf::F, Ks::F
     node.Ks = Param(Ks, bounds=node.Ks.bounds::Tuple)
 end
 
-
-# function update_state!(node::HyModNode, Sm, Sf1, Sf2, Sf3, Ss1, Q)::Nothing
-#     append!(node.Sm, Sm)
-#     append!(node.Sf1, Sf1)
-#     append!(node.Sf2, Sf2)
-#     append!(node.Sf3, Sf3)
-#     append!(node.Ss1, Ss1)
-#     append!(node.outflow, Q)
-
-#     return nothing
-# end
 function update_state!(
     node::HyModNode, ts::Int, Sm::F, Sf1::F, Sf2::F, Sf3::F, Ss1::F, Q::F
 )::Nothing where {F<:Float64}
diff --git a/src/Nodes/IHACRES/IHACRESNode.jl b/src/Nodes/IHACRES/IHACRESNode.jl
index 3c816df5..58a6a542 100644
--- a/src/Nodes/IHACRES/IHACRESNode.jl
+++ b/src/Nodes/IHACRES/IHACRESNode.jl
@@ -8,7 +8,7 @@ include("./components/flow.jl")
 abstract type IHACRESNode <: NetworkNode end
 
 
-Base.@kwdef mutable struct IHACRESBilinearNode{P, A<:AbstractFloat} <: IHACRESNode
+Base.@kwdef mutable struct IHACRESBilinearNode{P,A<:AbstractFloat} <: IHACRESNode
     const name::String
     const area::A
 
@@ -21,7 +21,7 @@ Base.@kwdef mutable struct IHACRESBilinearNode{P, A<:AbstractFloat} <: IHACRESNo
     b::P = Param(0.1, bounds=(1e-3, 0.1))  # slowflow storage coefficent == (1/tau_s)
 
     storage_coef::P = Param(2.9, bounds=(1e-10, 10.0))
-    alpha::P = Param(0.1, bounds=(1e-5, 1 - 1/10^9))
+    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
@@ -50,13 +50,13 @@ end
 
 
 function prep_state!(node::T, timesteps::Int64)::Nothing where {T<:IHACRESNode}
-    resize!(node.storage, timesteps+1)
+    resize!(node.storage, timesteps + 1)
     node.storage[2:end] .= 0.0
 
-    resize!(node.quick_store, timesteps+1)
+    resize!(node.quick_store, timesteps + 1)
     node.quick_store[2:end] .= 0.0
 
-    resize!(node.slow_store, timesteps+1)
+    resize!(node.slow_store, timesteps + 1)
     node.slow_store[2:end] .= 0.0
 
     node.outflow = fill(0.0, timesteps)
@@ -65,7 +65,7 @@ function prep_state!(node::T, timesteps::Int64)::Nothing where {T<:IHACRESNode}
     node.inflow = fill(0.0, timesteps)
     # node.level = fill(0.0, timesteps)
 
-    resize!(node.gw_store, timesteps+1)
+    resize!(node.gw_store, timesteps + 1)
     node.gw_store[2:end] .= 0.0
 
     return nothing
@@ -126,8 +126,8 @@ end
 Create a IHACRES node that adopts the bilinear CMD module.
 """
 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)
+    a::Float64, b::Float64, s_coef::Float64, alpha::Float64,
+    store::Float64, quick::Float64, slow::Float64, gw_store::Float64)
     n = create_node(IHACRESBilinearNode, name, area)
 
     update_params!(n, d, d2, e, f, a, b, s_coef, alpha)
@@ -192,9 +192,9 @@ end
 
 """
     run_timestep!(
-        node::IHACRESBilinearNode, climate::Climate, timestep::Int,
-        inflow::Float64, extraction::Float64, exchange::Float64
-    )
+        node::IHACRESNode, climate::Climate, timestep::Int64;
+        inflow=nothing, extraction=nothing, exchange=nothing
+    )::Float64
 
 Run the given IHACRESBilinearNode for a timestep.
 """
@@ -310,12 +310,12 @@ end
 
 
 function run_node_with_temp!(sn::StreamfallNetwork, nid::Int64, climate::Climate;
-                             inflow=nothing, extraction=nothing, exchange=nothing)
+    inflow=nothing, extraction=nothing, exchange=nothing)
     node = sn[nid]
     timesteps = sim_length(climate)
     @inbounds for ts in 1:timesteps
         run_node_with_temp!(node, climate, ts;
-                            inflow=inflow, extraction=extraction, exchange=exchange)
+            inflow=inflow, extraction=extraction, exchange=exchange)
     end
 
     return node.outflow, node.level
@@ -323,11 +323,11 @@ end
 
 
 function run_node_with_temp!(node::IHACRESBilinearNode, climate::Climate;
-                             inflow=nothing, extraction=nothing, exchange=nothing)
+    inflow=nothing, extraction=nothing, exchange=nothing)
     timesteps = sim_length(climate)
     @inbounds for ts in 1:timesteps
         run_node_with_temp!(node, climate, ts;
-                            inflow=inflow, extraction=extraction, exchange=exchange)
+            inflow=inflow, extraction=extraction, exchange=exchange)
     end
 
     return node.outflow, node.level
@@ -335,7 +335,7 @@ end
 
 
 function run_node_with_temp!(node::IHACRESBilinearNode, climate::Climate, timestep::Int64;
-                             inflow=nothing, extraction=nothing, exchange=nothing)
+    inflow=nothing, extraction=nothing, exchange=nothing)
     ts = timestep
     P, T = climate_values(node, climate, ts)
     i = timestep_value(ts, node.name, "_inflow", inflow)
@@ -346,36 +346,39 @@ end
 
 
 """
-    run_node_with_temp!(s_node::IHACRESBilinearNode,
-                        rain::Float64,
-                        temp::Float64,
-                        inflow::Float64,
-                        ext::Float64,
-                        gw_exchange::Float64=0.0;
-                        current_store=nothing,
-                        quick_store=nothing,
-                        slow_store=nothing
-                        gw_store=nothing)::Tuple{Float64, Float64}
+    run_node_with_temp!(
+        s_node::IHACRESBilinearNode,
+        rain::Float64,
+        temp::Float64,
+        inflow::Float64,
+        ext::Float64,
+        gw_exchange::Float64=0.0;
+        current_store=nothing,
+        quick_store=nothing,
+        slow_store=nothing,
+        gw_store=nothing
+    )::Tuple{Float64,Float64}
 
 Run node with temperature data to calculate outflow and update state.
 """
-function run_node_with_temp!(s_node::IHACRESBilinearNode,
-                             rain::Float64,
-                             temp::Float64,
-                             inflow::Float64,
-                             ext::Float64,
-                             gw_exchange::Float64=0.0;
-                             current_store=nothing,
-                             quick_store=nothing,
-                             slow_store=nothing,
-                             gw_store=nothing)::Tuple{Float64, Float64}
+function run_node_with_temp!(
+    s_node::IHACRESBilinearNode,
+    rain::Float64,
+    temp::Float64,
+    inflow::Float64,
+    ext::Float64,
+    gw_exchange::Float64=0.0;
+    current_store=nothing,
+    quick_store=nothing,
+    slow_store=nothing,
+    gw_store=nothing
+)::Tuple{Float64,Float64}
     current_store = s_node.storage[end]
     quick_store = s_node.quick_store[end]
     slow_store = s_node.slow_store[end]
     gw_store = s_node.gw_store[end]
 
     (mf, e_rainfall, recharge) = calc_ft_interim_cmd(
-        interim_results,
         current_store,
         rain,
         s_node.d,
diff --git a/src/Nodes/IHACRES/components/cmd.jl b/src/Nodes/IHACRES/components/cmd.jl
index c21597af..2975e48c 100644
--- a/src/Nodes/IHACRES/components/cmd.jl
+++ b/src/Nodes/IHACRES/components/cmd.jl
@@ -1,5 +1,8 @@
 """
-    calc_cmd(prev_cmd::Float64, rainfall::Float64, et::Float64, effective_rainfall::Float64, recharge::Float64)::Float64
+    calc_cmd(
+        prev_cmd::Float64, rainfall::Float64, et::Float64,
+        effective_rainfall::Float64, recharge::Float64
+    )::Float64
 
 Calculate Catchment Moisture Deficit (CMD).
 
@@ -75,7 +78,7 @@ function calc_trig_interim_cmd(cmd::Float64, param_d::Float64, rainfall::Float64
 end
 
 """
-    calc_ft_interim_cmd(cmd::Float64, rain::Float64, d::Float64, d2::Float64, alpha::Float64)
+    calc_ft_interim_cmd(cmd::Float64, rain::Float64, d::Float64, d2::Float64, alpha::Float64)::Tuple{Float64,Float64,Float64}
 
 Direct port of original Fortran implementation to calculate interim CMD (M_f).
 Calculates estimates of effective rainfall and recharge as a by-product.
@@ -90,9 +93,9 @@ Calculates estimates of effective rainfall and recharge as a by-product.
 # 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}
+function calc_ft_interim_cmd(cmd::Float64, rain::Float64, d::Float64, d2::Float64, alpha::Float64)::Tuple{Float64,Float64,Float64}
     # Initialize variables
-    d2 = d * d2_factor
+    d2 = d * d2
     tmp_cmd = cmd
     e_rain = 0.0
     recharge = 0.0
diff --git a/src/Nodes/IHACRES/components/flow.jl b/src/Nodes/IHACRES/components/flow.jl
index d736721e..42a7e245 100644
--- a/src/Nodes/IHACRES/components/flow.jl
+++ b/src/Nodes/IHACRES/components/flow.jl
@@ -60,8 +60,8 @@ Tuple of (quickflow, slowflow, outflow) in ML/day
     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}
+    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)
@@ -77,19 +77,19 @@ end
 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
+- `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}
+    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
@@ -141,8 +141,8 @@ Unit Hydrograph module ported from Fortran.
 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}
+    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
@@ -190,9 +190,9 @@ function calc_ft_level(outflow::Float64, level_params::Vector{Float64})::Float64
     end
 
     level = (exp(p1) * outflow^p2) *
-        (1.0 / (
-            1.0 + ((outflow / p3)^p4)^(p5/p4) * exp(p6 / (1.0 + exp(-p7*p8)) * outflow^p7)
-        ))
+            (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)
 
diff --git a/src/Nodes/IHACRES/components/p_and_et.jl b/src/Nodes/IHACRES/components/p_and_et.jl
index 25526a29..2c6b02b0 100644
--- a/src/Nodes/IHACRES/components/p_and_et.jl
+++ b/src/Nodes/IHACRES/components/p_and_et.jl
@@ -1,19 +1,22 @@
 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)
+    calc_effective_rainfall(
+        rainfall::Float64, cmd::Float64, d::Float64,
+        d2::Float64; n::Float64=0.1
+    )::Float64
 
 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.
+  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]"
+  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
 
@@ -27,8 +30,10 @@ Estimate effective rainfall.
 # Returns
 - Effective rainfall value
 """
-function calc_effective_rainfall(rainfall::Float64, cmd::Float64, d::Float64,
-                               d2::Float64, n::Float64=0.1)::Float64
+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
@@ -43,7 +48,7 @@ function calc_effective_rainfall(rainfall::Float64, cmd::Float64, d::Float64,
 end
 
 """
-    calc_ET_from_E(e::Float64, evap::Float64, Mf::Float64, f::Float64, d::Float64)
+    calc_ET_from_E(e::F, evap::F, Mf::F, f::F, d::F)::F where {F<:Float64}
 
 Calculate evapotranspiration from evaporation.
 
@@ -57,30 +62,19 @@ Calculate evapotranspiration from evaporation.
 # Returns
 - Estimate of ET
 """
-function calc_ET_from_E(e::Float64, evap::Float64, Mf::Float64, f::Float64,
-                       d::Float64)::Float64
+function calc_ET_from_E(e::F, evap::F, Mf::F, f::F, d::F)::F where {F<:Float64}
     param_g = f * d
     et = e * evap
 
     if Mf > param_g
-        et *= min(1.0, exp((1.0 - Mf/param_g)*2.0))
+        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)
+    calc_ET_from_T(e::F, T::F, Mf::F, f::F, d::F)::F where {F<:Float64}
 
 Calculate evapotranspiration based on temperature data.
 
@@ -96,10 +90,9 @@ water availability for plant transpiration.
 - `d`: flow threshold factor
 
 # Returns
-- Estimate of ET from temperature (for catchment area)
+Estimate of ET from temperature (for catchment area)
 """
-function calc_ET_from_T(e::Float64, T::Float64, Mf::Float64, f::Float64,
-                       d::Float64)::Float64
+function calc_ET_from_T(e::F, T::F, Mf::F, f::F, d::F)::F where {F<:Float64}
     # temperature can be negative, so we have a min cap of 0.0
     if T <= 0.0
         return 0.0
diff --git a/src/metrics.jl b/src/metrics.jl
index c10e6e69..dad14ec0 100644
--- a/src/metrics.jl
+++ b/src/metrics.jl
@@ -110,7 +110,7 @@ end
 
 
 """
-Applies split meta metric approach
+Applies split meta metric approach.
 
 If using with other macros such as `@normalize` or `@bound`,
 these must come first.
@@ -136,25 +136,45 @@ macro split(metric, obs, sim, n, agg_func=mean)
 end
 
 
-"""The Nash-Sutcliffe Efficiency score"""
-NSE(obs, sim) = 1.0 - sum((obs .- sim) .^ 2) / sum((obs .- mean(obs)) .^ 2)
+"""
+    NSE(obs, sim)
 
+The Nash-Sutcliffe Efficiency score
+"""
+function NSE(obs, sim)
+    return 1.0 - sum((obs .- sim) .^ 2) / sum((obs .- mean(obs)) .^ 2)
+end
 
-"""Normalized Nash-Sutcliffe Efficiency score (bounded between 0 and 1).
+
+"""
+    NNSE(obs, sim)
+
+Normalized Nash-Sutcliffe Efficiency score (bounded between 0 and 1).
 
 # References
 1. Nossent, J., Bauwens, W., 2012.
     Application of a normalized Nash-Sutcliffe efficiency to improve the accuracy of the Sobol’ sensitivity analysis of a hydrological model.
     EGU General Assembly Conference Abstracts 237.
 """
-NNSE(obs, sim) = 1.0 / (2.0 - NSE(obs, sim))
+function NNSE(obs, sim)
+    return 1.0 / (2.0 - NSE(obs, sim))
+end
 
 
-"""Root Mean Square Error"""
-RMSE(obs, sim) = (sum((sim .- obs) .^ 2) / length(sim))^0.5
+"""
+    RMSE(obs, sim)
+
+Root Mean Square Error
+"""
+function RMSE(obs, sim)
+    return (sum((sim .- obs) .^ 2) / length(sim))^0.5
+end
+
 
+"""
+    R2(obs, sim)::Float64
 
-"""Coefficient of determination (R^2)
+Coefficient of determination (R²)
 
 Aliases `NSE()`
 """
@@ -163,12 +183,15 @@ function R2(obs, sim)::Float64
 end
 
 
-"""Adjusted R^2
+"""
+    ADJ_R2(obs, sim, p::Int64)::Float64
+
+Adjusted R²
 
 # Arguments
 - `obs::Vector` : observations
 - `sim::Vector` : modeled results
-- `p::Int` : number of explanatory variables
+- `p::Int64` : number of explanatory variables
 """
 function ADJ_R2(obs, sim, p::Int64)::Float64
     n = length(obs)
@@ -179,12 +202,16 @@ end
 
 
 """
+    MAE(obs, sim)
+
 Mean Absolute Error
 """
 MAE(obs, sim) = mean(abs.(sim .- obs))
 
 
 """
+    ME(obs, sim)
+
 Mean Error
 """
 ME(obs, sim) = mean(sim .- obs)
@@ -249,6 +276,8 @@ end
 
 
 """
+    relative_skill_score(Sb::Float64, Sm::Float64)::Float64
+
 Relative Skill Score.
 
 Provides an indication of model performance relative to a known benchmark score.
@@ -266,13 +295,13 @@ Suitable for use with least-squares approaches that provide skill scores ranging
     Catchment hydrology/Modelling approaches.
     https://doi.org/10.5194/hess-2019-327
 """
-function relative_skill_score(Sb::Float64, Sm::Float64)
+function relative_skill_score(Sb::Float64, Sm::Float64)::Float64
     return (Sm - Sb) / (1.0 - Sb)
 end
 
 
 """
-    NSE_logbias(obs, sim; metric::Function=NSE, bias_threshold::Float64=5.0, shape::Float64=2.5)
+    NSE_logbias(obs, sim; metric=NSE, bias_threshold=5.0, shape=2.5)
 
 The NSE_logbias meta-metric provides a weighted combination of a least-squares approach and a
 logarithmic function of bias. The metric penalizes predictions with an overall bias above a
@@ -313,7 +342,7 @@ end
 
 
 """
-    KGE(obs::Vector, sim::Vector; scaling::Tuple=nothing)::Float64
+    KGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
 
 Calculate the 2009 Kling-Gupta Efficiency (KGE) metric.
 
@@ -348,11 +377,7 @@ Note: Although similar, NSE and KGE cannot be directly compared.
     Hydrology and Earth System Sciences 23, 2601–2614.
     https://doi.org/10.5194/hess-23-2601-2019
 """
-function KGE(obs, sim; scaling=nothing)::Float64
-    if isnothing(scaling)
-        scaling = (1, 1, 1)
-    end
-
+function KGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
     # Correlation
     r = Statistics.cor(obs, sim)
     if isnan(r)
@@ -375,7 +400,10 @@ function KGE(obs, sim; scaling=nothing)::Float64
 end
 
 
-"""Bounded KGE, bounded between -1 and 1.
+"""
+    BKGE(obs, sim)::Float64
+
+Bounded KGE, bounded between -1 and 1.
 
 # Arguments
 - `obs::Vector` : observations
@@ -387,18 +415,24 @@ function BKGE(obs, sim)::Float64
 end
 
 
-"""Normalized KGE between 0 and 1.
+"""
+    NKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
+
+Normalized KGE between 0 and 1.
 
 # Arguments
 - `obs::Vector` : observations
 - `sim::Vector` : modeled results
 """
-function NKGE(obs, sim; scaling=nothing)::Float64
+function NKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
     return 1 / (2 - KGE(obs, sim; scaling=scaling))
 end
 
 
-"""Calculate the modified KGE metric (2012).
+"""
+    mKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
+
+Calculate the modified KGE metric (2012).
 
 Also known as KGE prime (KGE').
 
@@ -433,11 +467,7 @@ This is to:
     Hydrology and Earth System Sciences 22, 4583–4591.
     https://doi.org/10.5194/hess-22-4583-2018
 """
-function mKGE(obs, sim; scaling=nothing)::Float64
-    if isnothing(scaling)
-        scaling = (1, 1, 1)
-    end
-
+function mKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
     # Timing (Pearson's correlation)
     r = Statistics.cor(obs, sim)
     if isnan(r)
@@ -475,9 +505,7 @@ function mKGE(obs, sim; scaling=nothing)::Float64
         β = μ_s / μ_o
     end
 
-    rs = scaling[1]
-    βs = scaling[2]
-    γs = scaling[3]
+    rs, βs, γs = scaling
 
     mod_kge = 1.0 - sqrt(rs * (r - 1)^2 + βs * (β - 1)^2 + γs * (γ - 1)^2)
 
@@ -485,34 +513,41 @@ function mKGE(obs, sim; scaling=nothing)::Float64
 end
 
 
-"""Bounded modified KGE between -1 and 1.
+"""
+    BmKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
+
+Bounded modified KGE between -1 and 1.
 
 # Arguments
 - `obs::Vector` : observations
 - `sim::Vector` : modeled results
 """
-function BmKGE(obs, sim; scaling=nothing)::Float64
+function BmKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
     mkge = mKGE(obs, sim; scaling=scaling)
     return mkge / (2 - mkge)
 end
 
 
-"""Normalized modified KGE between 0 and 1.
+"""
+    NmKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
+
+Normalized modified KGE between 0 and 1.
 
 # Arguments
 - `obs::Vector` : observations
 - `sim::Vector` : modeled results
 """
-function NmKGE(obs, sim; scaling=nothing)::Float64
+function NmKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
     return 1 / (2 - mKGE(obs, sim; scaling=scaling))
 end
 
 
 """
+    mean_NmKGE(obs, sim; scaling=(1.0, 1.0, 1.0), ϵ=1e-2)
+
 Mean Inverse NmKGE
 
-Said to produce better fits for low-flow indices
-compared to mKGE (see [1]).
+Said to produce better fits for low-flow indices compared to mKGE (see [1]).
 
 # Arguments
 - `obs::Vector` : observations
@@ -527,15 +562,23 @@ compared to mKGE (see [1]).
     Hydrological Sciences Journal 62, 1149–1166.
     https://doi.org/10.1080/02626667.2017.1308511
 """
-mean_NmKGE(obs, sim; scaling=nothing, ϵ=1e-2) = mean([Streamfall.NmKGE(obs, sim; scaling=scaling), Streamfall.NmKGE(1.0 ./ (obs .+ ϵ), 1.0 ./ (sim .+ ϵ); scaling=scaling)])
+function mean_NmKGE(obs, sim; scaling=(1.0, 1.0, 1.0), ϵ=1e-2)
+    return mean([
+        Streamfall.NmKGE(obs, sim; scaling=scaling),
+        Streamfall.NmKGE(1.0 ./ (obs .+ ϵ), 1.0 ./ (sim .+ ϵ); scaling=scaling)
+    ])
+end
 
 
-"""Calculate the non-parametric Kling-Gupta Efficiency (KGE) metric.
+"""
+    npKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
+
+Calculate the non-parametric Kling-Gupta Efficiency (KGE) metric.
 
 # Arguments
-- `obs::Vector` : observations
-- `sim::Vector` : modeled
-- `scaling::Tuple` : scaling factors for timing (s), variability (α), magnitude (β)
+- `obs` : observations
+- `sim` : modeled
+- `scaling` : scaling factors for timing (s), variability (α), magnitude (β)
 
 # References
 1. Pool, S., Vis, M., Seibert, J., 2018.
@@ -544,12 +587,8 @@ mean_NmKGE(obs, sim; scaling=nothing, ϵ=1e-2) = mean([Streamfall.NmKGE(obs, sim
     https://doi.org/10.1080/02626667.2018.1552002
 
 """
-function npKGE(obs, sim; scaling=nothing)::Float64
-    if isnothing(scaling)
-        scaling = (1, 1, 1)
-    end
-
-    # flow duration curves
+function npKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
+    # Flow duration curves
     μ_s = mean(sim)
     if μ_s == 0.0
         fdc_sim = repeat([0.0], length(sim))
@@ -581,9 +620,7 @@ function npKGE(obs, sim; scaling=nothing)::Float64
         r = 1.0  # can occur if identical sequences are used (e.g., 0 flows)
     end
 
-    rs = scaling[1]
-    αs = scaling[2]
-    βs = scaling[3]
+    rs, αs, βs = scaling
 
     kge = 1 - sqrt(rs * (r - 1)^2 + αs * (α - 1)^2 + βs * (β - 1)^2)
 
@@ -591,13 +628,16 @@ function npKGE(obs, sim; scaling=nothing)::Float64
 end
 
 
-"""Bounded non-parametric KGE between -1 and 1.
+"""
+    BnpKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
+
+Bounded non-parametric KGE between -1 and 1.
 
 # Arguments
 - `obs::Vector` : observations
 - `sim::Vector` : modeled results
 """
-function BnpKGE(obs, sim; scaling=nothing)::Float64
+function BnpKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
     npkge = npKGE(obs, sim; scaling=scaling)
     return npkge / (2 - npkge)
 end
@@ -606,18 +646,20 @@ end
 """Normalized non-parametric KGE between 0 and 1.
 
 # Arguments
-- `obs::Vector` : observations
-- `sim::Vector` : modeled results
+- `obs` : observations
+- `sim` : modeled results
 """
-function NnpKGE(obs, sim; scaling=nothing)::Float64
+function NnpKGE(obs, sim; scaling=(1.0, 1.0, 1.0))::Float64
     return 1 / (2 - npKGE(obs, sim; scaling=scaling))
 end
 
 
-"""Liu Mean Efficiency metric (LME).
+"""
+    LME(obs, sim)::Float64
+
+Liu Mean Efficiency metric (LME).
 
-Reformulation of the KGE metric said to be advantageous for capturing extreme
-flow events.
+Reformulation of the KGE metric said to be advantageous for capturing extreme flow events.
 
 # Arguments
 - `obs::Vector` : observations
@@ -659,7 +701,10 @@ function naive_split_metric(obs::Vector, sim::Vector, n_members::Int, metric::Fu
 end
 
 
-"""Naive approach to split metrics.
+"""
+    naive_split_metric(obs, sim; n_members::Int=365, metric::Function=NNSE, comb_method::Function=mean)
+
+Naive approach to split metrics.
 
 Split metrics are a meta-objective optimization approach which "splits" data
 into subperiods. The objective function is calculated for each subperiod and
@@ -690,20 +735,19 @@ end
 
 
 """
-    inverse_metric(obs, sim, metric; comb_method::Function=mean, ϵ=1e-2)
+    inverse_metric(obs, sim, metric::Function; comb_method::Function=mean, ϵ=1e-2)
 
-A meta-objective function which combines the performance of the
-given metric as applied to the discharge and the inverse of the
-discharge.
+A meta-objective function which combines the performance of the given metric as applied to
+the discharge and the inverse of the discharge.
 
 By default, the combination method is to take the mean.
 
 # Arguments
-- obs : observed
-- sim : modeled results
-- `metric::Function` : objective function
-- `comb_method::Function` : mean
-- ϵ : offset value to use (enables use with zero-flow time steps), defaults to 1e-2
+- `obs` : observed
+- `sim` : modeled results
+- `metric` : objective function
+- `comb_method` : Method to combine outputs (default: `mean`)
+- `ϵ` : offset value to use (enables use with zero-flow time steps), defaults to 1e-2
 
 # References
 1. Garcia, F., Folton, N., Oudin, L., 2017.
@@ -718,6 +762,8 @@ end
 
 
 """
+    skill_score(model_score, benchmark_score)
+
 Allows comparison of any model compared against a pre-defined benchmark, assuming both scores were obtained with
 the same objective function.