diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl
index dc9613688..23b83d8ca 100644
--- a/core/src/allocation_optim.jl
+++ b/core/src/allocation_optim.jl
@@ -311,7 +311,7 @@ function get_basin_data(
     (; mean_input_flows) = allocation
     @assert node_id.type == NodeType.Basin
     influx = mean_input_flows[(node_id, node_id)][]
-    storage_basin = basin.current_storage[parent(u)][node_id.idx]
+    storage_basin = basin.current_properties.current_storage[parent(u)][node_id.idx]
     control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control)
     if isempty(control_inneighbors)
         level_demand_idx = 0
diff --git a/core/src/bmi.jl b/core/src/bmi.jl
index 8e5c9960c..c02b27278 100644
--- a/core/src/bmi.jl
+++ b/core/src/bmi.jl
@@ -33,9 +33,9 @@ end
 function BMI.get_value_ptr(model::Model, name::AbstractString)::AbstractVector{Float64}
     (; u, p) = model.integrator
     if name == "basin.storage"
-        p.basin.current_storage[parent(u)]
+        p.basin.current_properties.current_storage[parent(u)]
     elseif name == "basin.level"
-        p.basin.current_level[parent(u)]
+        p.basin.current_properties.current_level[parent(u)]
     elseif name == "basin.infiltration"
         p.basin.vertical_flux.infiltration
     elseif name == "basin.drainage"
diff --git a/core/src/callback.jl b/core/src/callback.jl
index b1c2b3ad3..ade71dc3f 100644
--- a/core/src/callback.jl
+++ b/core/src/callback.jl
@@ -316,8 +316,10 @@ function update_cumulative_flows!(u, t, integrator)::Nothing
     end
 
     # Update the Basin concentrations again based on the removed mass
-    basin.concentration_state .= basin.mass ./ basin.current_storage[parent(u)]
-    basin.storage_prev .= basin.current_storage[parent(u)]
+    basin.concentration_state .=
+        basin.mass ./ basin.current_properties.current_storage[parent(u)]
+    basin.storage_prev .= basin.current_properties.current_storage[parent(u)]
+    basin.level_prev .= basin.current_properties.current_level[parent(u)]
 
     return nothing
 end
@@ -361,8 +363,8 @@ function save_basin_state(u, t, integrator)
     (; p) = integrator
     (; basin) = p
     du = get_du(integrator)
-    current_storage = basin.current_storage[parent(du)]
-    current_level = basin.current_level[parent(du)]
+    current_storage = basin.current_properties.current_storage[parent(du)]
+    current_level = basin.current_properties.current_level[parent(du)]
     water_balance!(du, u, p, t)
     SavedBasinState(; storage = copy(current_storage), level = copy(current_level), t)
 end
@@ -448,7 +450,7 @@ function check_water_balance_error!(
     (; u, p, t) = integrator
     (; basin, water_balance_abstol, water_balance_reltol) = p
     errors = false
-    current_storage = basin.current_storage[parent(u)]
+    current_storage = basin.current_properties.current_storage[parent(u)]
 
     # The initial storage is irrelevant for the storage rate and can only cause
     # floating point truncation errors
@@ -514,7 +516,8 @@ end
 
 function check_negative_storage(u, t, integrator)::Nothing
     (; basin) = integrator.p
-    (; node_id, current_storage) = basin
+    (; node_id, current_properties) = basin
+    (; current_storage) = current_properties
     du = get_du(integrator)
     set_current_basin_properties!(du, u, integrator.p, t)
     current_storage = current_storage[parent(du)]
@@ -737,7 +740,7 @@ end
 function update_subgrid_level!(integrator)::Nothing
     (; p) = integrator
     du = get_du(integrator)
-    basin_level = p.basin.current_level[parent(du)]
+    basin_level = p.basin.current_properties.current_level[parent(du)]
     subgrid = integrator.p.subgrid
     for (i, (index, interp)) in enumerate(zip(subgrid.basin_index, subgrid.interpolations))
         subgrid.level[i] = interp(basin_level[index])
@@ -823,7 +826,7 @@ update_userd_conc!(integrator)::Nothing =
 function update_allocation!(integrator)::Nothing
     (; p, t, u) = integrator
     (; allocation, basin) = p
-    (; current_storage) = basin
+    (; current_storage) = basin.current_properties
     (; allocation_models, mean_input_flows, mean_realized_flows) = allocation
 
     # Make sure current storages are up to date
diff --git a/core/src/model.jl b/core/src/model.jl
index d4aa9af1e..92d7ac3a2 100644
--- a/core/src/model.jl
+++ b/core/src/model.jl
@@ -119,6 +119,11 @@ function Model(config::Config)::Model
     # Synchronize level with storage
     set_current_basin_properties!(du0, u0, parameters, t0)
 
+    # Previous level is used to estimate the minimum level that was attained during a time step
+    # in limit_flow!
+    parameters.basin.level_prev .=
+        parameters.basin.current_properties.current_level[parent(u0)]
+
     saveat = convert_saveat(config.solver.saveat, t_end)
     saveat isa Float64 && push!(tstops, range(0, t_end; step = saveat))
     tstops = sort(unique(vcat(tstops...)))
diff --git a/core/src/parameter.jl b/core/src/parameter.jl
index 1f9a8f62d..97bc1cab3 100644
--- a/core/src/parameter.jl
+++ b/core/src/parameter.jl
@@ -1,3 +1,10 @@
+
+# Universal reduction factor threshold for the low storage factor
+const LOW_STORAGE_THRESHOLD = 10.0
+
+# Universal reduction factor threshold for the minimum upstream level of UserDemand nodes
+const USER_DEMAND_MIN_LEVEL_THRESHOLD = 0.1
+
 const SolverStats = @NamedTuple{
     time::Float64,
     rhs_calls::Int,
@@ -300,6 +307,23 @@ In-memory storage of saved instantaneous storages and levels for writing to resu
     t::Float64
 end
 
+"""
+Caches of current basin properties
+"""
+struct CurrentBasinProperties
+    current_storage::Cache
+    # Low storage factor for reducing flows out of drying basins
+    # given the current storages
+    current_low_storage_factor::Cache
+    current_level::Cache
+    current_area::Cache
+    current_cumulative_precipitation::Cache
+    current_cumulative_drainage::Cache
+    function CurrentBasinProperties(n)
+        new((cache(n) for _ in 1:6)...)
+    end
+end
+
 """
 Requirements:
 
@@ -332,11 +356,7 @@ end
     cumulative_precipitation_saveat::Vector{Float64} = zeros(length(node_id))
     cumulative_drainage_saveat::Vector{Float64} = zeros(length(node_id))
     # Cache this to avoid recomputation
-    current_storage::Cache = cache(length(node_id))
-    current_level::Cache = cache(length(node_id))
-    current_area::Cache = cache(length(node_id))
-    current_cumulative_precipitation::Cache = cache(length(node_id))
-    current_cumulative_drainage::Cache = cache(length(node_id))
+    current_properties::CurrentBasinProperties = CurrentBasinProperties(length(node_id))
     # Discrete values for interpolation
     storage_to_level::Vector{
         LinearInterpolationIntInv{
@@ -355,6 +375,8 @@ end
     # Data source for concentration updates
     concentration_time::StructVector{BasinConcentrationV1, D, Int}
 
+    # Level for each Basin at the previous time step
+    level_prev::Vector{Float64} = zeros(length(node_id))
     # Concentrations
     # Config setting to enable/disable evaporation of mass
     evaporate_mass::Bool = true
diff --git a/core/src/read.jl b/core/src/read.jl
index 2a9ae9bf5..749554c3b 100644
--- a/core/src/read.jl
+++ b/core/src/read.jl
@@ -574,13 +574,10 @@ end
 function Basin(db::DB, config::Config, graph::MetaGraph)::Basin
     node_id = get_ids(db, "Basin")
     n = length(node_id)
-    current_level = cache(n)
-    current_area = cache(n)
 
     evaporate_mass = config.solver.evaporate_mass
     precipitation = zeros(n)
     potential_evaporation = zeros(n)
-    evaporation = zeros(n)
     drainage = zeros(n)
     infiltration = zeros(n)
     table = (; precipitation, potential_evaporation, drainage, infiltration)
@@ -686,8 +683,6 @@ function Basin(db::DB, config::Config, graph::MetaGraph)::Basin
         inflow_ids = [collect(inflow_ids(graph, id)) for id in node_id],
         outflow_ids = [collect(outflow_ids(graph, id)) for id in node_id],
         vertical_flux,
-        current_level,
-        current_area,
         storage_to_level,
         level_to_area,
         demand,
diff --git a/core/src/solve.jl b/core/src/solve.jl
index 3863dc633..6d306fdf3 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -8,14 +8,17 @@ function water_balance!(
     t::Number,
 )::Nothing
     (; basin, pid_control) = p
+    (; current_storage, current_low_storage_factor, current_level) =
+        basin.current_properties
 
     du .= 0.0
 
     # Ensures current_* vectors are current
     set_current_basin_properties!(du, u, p, t)
 
-    current_storage = basin.current_storage[parent(du)]
-    current_level = basin.current_level[parent(du)]
+    current_storage = current_storage[parent(du)]
+    current_low_storage_factor = current_low_storage_factor[parent(du)]
+    current_level = current_level[parent(du)]
 
     # Notes on the ordering of these formulations:
     # - Continuous control can depend on flows (which are not continuously controlled themselves),
@@ -27,7 +30,7 @@ function water_balance!(
     update_vertical_flux!(basin, du)
 
     # Formulate intermediate flows (non continuously controlled)
-    formulate_flows!(du, p, t, current_storage, current_level)
+    formulate_flows!(du, p, t, current_storage, current_low_storage_factor, current_level)
 
     # Compute continuous control
     formulate_continuous_control!(du, p, t)
@@ -38,6 +41,7 @@ function water_balance!(
         p,
         t,
         current_storage,
+        current_low_storage_factor,
         current_level;
         continuous_control_type = ContinuousControlType.Continuous,
     )
@@ -51,6 +55,7 @@ function water_balance!(
         p,
         t,
         current_storage,
+        current_low_storage_factor,
         current_level;
         continuous_control_type = ContinuousControlType.PID,
     )
@@ -79,17 +84,19 @@ function set_current_basin_properties!(
     t::Number,
 )::Nothing
     (; basin) = p
+    (; current_properties, cumulative_precipitation, cumulative_drainage, vertical_flux) =
+        basin
     (;
         current_storage,
+        current_low_storage_factor,
         current_level,
         current_area,
         current_cumulative_precipitation,
         current_cumulative_drainage,
-        cumulative_precipitation,
-        cumulative_drainage,
-        vertical_flux,
-    ) = basin
+    ) = current_properties
+
     current_storage = current_storage[parent(du)]
+    current_low_storage_factor = current_low_storage_factor[parent(du)]
     current_level = current_level[parent(du)]
     current_area = current_area[parent(du)]
     current_cumulative_precipitation = current_cumulative_precipitation[parent(du)]
@@ -107,7 +114,9 @@ function set_current_basin_properties!(
 
     formulate_storages!(current_storage, du, u, p, t)
 
-    for (i, s) in enumerate(current_storage)
+    for (id, s) in zip(basin.node_id, current_storage)
+        i = id.idx
+        current_low_storage_factor[i] = reduction_factor(s, LOW_STORAGE_THRESHOLD)
         current_level[i] = get_level_from_storage(basin, i, s)
         current_area[i] = basin.level_to_area[i](current_level[i])
     end
@@ -144,7 +153,8 @@ function formulate_storage!(
     basin::Basin,
     du::ComponentVector,
 )
-    (; current_cumulative_precipitation, current_cumulative_drainage) = basin
+    (; current_cumulative_precipitation, current_cumulative_drainage) =
+        basin.current_properties
 
     current_cumulative_precipitation = current_cumulative_precipitation[parent(du)]
     current_cumulative_drainage = current_cumulative_drainage[parent(du)]
@@ -185,7 +195,8 @@ Smoothly let the evaporation flux go to 0 when at small water depths
 Currently at less than 0.1 m.
 """
 function update_vertical_flux!(basin::Basin, du::AbstractVector)::Nothing
-    (; vertical_flux, current_level, current_area) = basin
+    (; vertical_flux, current_properties) = basin
+    (; current_level, current_area) = current_properties
     current_level = current_level[parent(du)]
     current_area = current_area[parent(du)]
 
@@ -211,7 +222,7 @@ function set_error!(pid_control::PidControl, p::Parameters, du::ComponentVector,
     (; basin) = p
     (; listen_node_id, target, error) = pid_control
     error = error[parent(du)]
-    current_level = basin.current_level[parent(du)]
+    current_level = basin.current_properties.current_level[parent(du)]
 
     for i in eachindex(listen_node_id)
         listened_node_id = listen_node_id[i]
@@ -229,7 +240,7 @@ function formulate_pid_control!(
 )::Nothing
     (; basin) = p
     (; node_id, active, target, listen_node_id, error) = pid_control
-    (; current_area) = basin
+    (; current_area) = basin.current_properties
 
     current_area = current_area[parent(du)]
     error = error[parent(du)]
@@ -315,7 +326,7 @@ function formulate_flow!(
     user_demand::UserDemand,
     p::Parameters,
     t::Number,
-    current_storage::Vector,
+    current_low_storage_factor::Vector,
     current_level::Vector,
 )::Nothing
     (; allocation) = p
@@ -348,14 +359,14 @@ function formulate_flow!(
 
         # Smoothly let abstraction go to 0 as the source basin dries out
         inflow_id = inflow_edge.edge[1]
-        factor_basin = low_storage_factor(current_storage, inflow_id, 10.0)
+        factor_basin = get_low_storage_factor(current_low_storage_factor, inflow_id)
         q *= factor_basin
 
         # Smoothly let abstraction go to 0 as the source basin
         # level reaches its minimum level
         source_level = get_level(p, inflow_id, t, current_level)
         Δsource_level = source_level - min_level
-        factor_level = reduction_factor(Δsource_level, 0.1)
+        factor_level = reduction_factor(Δsource_level, USER_DEMAND_MIN_LEVEL_THRESHOLD)
         q *= factor_level
         du.user_demand_inflow[id.idx] = q
         du.user_demand_outflow[id.idx] = q * return_factor(t)
@@ -371,7 +382,7 @@ function formulate_flow!(
     linear_resistance::LinearResistance,
     p::Parameters,
     t::Number,
-    current_storage::Vector,
+    current_low_storage_factor::Vector,
     current_level::Vector,
 )::Nothing
     all_nodes_active = p.all_nodes_active[]
@@ -389,11 +400,10 @@ function formulate_flow!(
             q_unlimited = (h_a - h_b) / resistance[id.idx]
             q = clamp(q_unlimited, -max_flow_rate[id.idx], max_flow_rate[id.idx])
             q *= low_storage_factor_resistance_node(
-                current_storage,
+                current_low_storage_factor,
                 q,
                 inflow_id,
                 outflow_id,
-                10.0,
             )
             du.linear_resistance[id.idx] = q
         end
@@ -409,10 +419,9 @@ function formulate_flow!(
     tabulated_rating_curve::TabulatedRatingCurve,
     p::Parameters,
     t::Number,
-    current_storage::Vector,
+    current_low_storage_factor::Vector,
     current_level::Vector,
 )::Nothing
-    (; basin) = p
     all_nodes_active = p.all_nodes_active[]
     (; node_id, active, table) = tabulated_rating_curve
 
@@ -428,7 +437,7 @@ function formulate_flow!(
         Δh = h_a - h_b
 
         if active[id.idx] || all_nodes_active
-            factor = low_storage_factor(current_storage, inflow_id, 10.0)
+            factor = get_low_storage_factor(current_low_storage_factor, inflow_id)
             q = factor * table[id.idx](h_a)
             q *= reduction_factor(Δh, 0.02)
             q *= reduction_factor(max_downstream_level - h_b, 0.02)
@@ -485,7 +494,7 @@ function formulate_flow!(
     manning_resistance::ManningResistance,
     p::Parameters,
     t::Number,
-    current_storage::Vector,
+    current_low_storage_factor::Vector,
     current_level::Vector,
 )::Nothing
     (;
@@ -540,11 +549,10 @@ function formulate_flow!(
 
         q = A / n * ∛(R_h^2) * relaxed_root(Δh / L, 1e-3)
         q *= low_storage_factor_resistance_node(
-            current_storage,
+            current_low_storage_factor,
             q,
             inflow_id,
             outflow_id,
-            10.0,
         )
         du.manning_resistance[id.idx] = q
     end
@@ -556,7 +564,7 @@ function formulate_flow!(
     pump::Pump,
     p::Parameters,
     t::Number,
-    current_storage::Vector,
+    current_low_storage_factor::Vector,
     current_level::Vector,
     continuous_control_type_::ContinuousControlType.T,
 )::Nothing
@@ -594,7 +602,7 @@ function formulate_flow!(
         src_level = get_level(p, inflow_id, t, current_level)
         dst_level = get_level(p, outflow_id, t, current_level)
 
-        factor = low_storage_factor(current_storage, inflow_id, 10.0)
+        factor = get_low_storage_factor(current_low_storage_factor, inflow_id)
         q = flow_rate * factor
 
         q *= reduction_factor(src_level - min_upstream_level, 0.02)
@@ -611,7 +619,7 @@ function formulate_flow!(
     outlet::Outlet,
     p::Parameters,
     t::Number,
-    current_storage::Vector,
+    current_low_storage_factor::Vector,
     current_level::Vector,
     continuous_control_type_::ContinuousControlType.T,
 )::Nothing
@@ -650,7 +658,7 @@ function formulate_flow!(
         dst_level = get_level(p, outflow_id, t, current_level)
 
         q = flow_rate
-        q *= low_storage_factor(current_storage, inflow_id, 10.0)
+        q *= get_low_storage_factor(current_low_storage_factor, inflow_id)
 
         # No flow of outlet if source level is lower than target level
         Δlevel = src_level - dst_level
@@ -669,6 +677,7 @@ function formulate_flows!(
     p::Parameters,
     t::Number,
     current_storage::Vector,
+    current_low_storage_factor::Vector,
     current_level::Vector;
     continuous_control_type::ContinuousControlType.T = ContinuousControlType.None,
 )::Nothing
@@ -681,22 +690,51 @@ function formulate_flows!(
         user_demand,
     ) = p
 
-    formulate_flow!(du, pump, p, t, current_storage, current_level, continuous_control_type)
+    formulate_flow!(
+        du,
+        pump,
+        p,
+        t,
+        current_low_storage_factor,
+        current_level,
+        continuous_control_type,
+    )
     formulate_flow!(
         du,
         outlet,
         p,
         t,
-        current_storage,
+        current_low_storage_factor,
         current_level,
         continuous_control_type,
     )
 
     if continuous_control_type == ContinuousControlType.None
-        formulate_flow!(du, linear_resistance, p, t, current_storage, current_level)
-        formulate_flow!(du, manning_resistance, p, t, current_storage, current_level)
-        formulate_flow!(du, tabulated_rating_curve, p, t, current_storage, current_level)
-        formulate_flow!(du, user_demand, p, t, current_storage, current_level)
+        formulate_flow!(
+            du,
+            linear_resistance,
+            p,
+            t,
+            current_low_storage_factor,
+            current_level,
+        )
+        formulate_flow!(
+            du,
+            manning_resistance,
+            p,
+            t,
+            current_low_storage_factor,
+            current_level,
+        )
+        formulate_flow!(
+            du,
+            tabulated_rating_curve,
+            p,
+            t,
+            current_low_storage_factor,
+            current_level,
+        )
+        formulate_flow!(du, user_demand, p, t, current_low_storage_factor, current_level)
     end
 end
 
@@ -721,6 +759,14 @@ function limit_flow!(
         allocation,
     ) = p
 
+    # The current storage and level based on the proposed u are used to estimate the lowest
+    # storage and level attained in the last time step to estimate whether there was an effect
+    # of reduction factors
+    du = get_du(integrator)
+    set_current_basin_properties!(du, u, p, t)
+    current_storage = basin.current_properties.current_storage[parent(u)]
+    current_level = basin.current_properties.current_level[parent(u)]
+
     # TabulatedRatingCurve flow is in [0, ∞) and can be inactive
     for (id, active) in zip(tabulated_rating_curve.node_id, tabulated_rating_curve.active)
         limit_flow!(
@@ -763,28 +809,60 @@ function limit_flow!(
         )
     end
 
-    # UserDemand inflow is in [0, total_demand] and can be inactive
-    for (id, active) in zip(user_demand.node_id, user_demand.active)
-        total_demand = sum(
-            get_demand(user_demand, id, priority_idx, t) for
-            priority_idx in eachindex(allocation.priorities)
-        )
+    # UserDemand inflow bounds depend on multiple aspects of the simulation
+    for (id, active, inflow_edge, demand_from_timeseries) in zip(
+        user_demand.node_id,
+        user_demand.active,
+        user_demand.inflow_edge,
+        user_demand.demand_from_timeseries,
+    )
+        min_flow_rate, max_flow_rate = if demand_from_timeseries
+            # Bounding the flow rate if the demand comes from a time series is hard
+            0, Inf
+        else
+            # The lower bound is estimated as the lowest inflow given the minimum values
+            # of the reduction factors involved (with a margin)
+            inflow_id = inflow_edge.edge[1]
+            factor_basin_min =
+                min_low_storage_factor(current_storage, basin.storage_prev, inflow_id)
+            factor_level_min = min_low_user_demand_level_factor(
+                current_level,
+                basin.level_prev,
+                user_demand.min_level,
+                id,
+                inflow_id,
+            )
+            allocated_total =
+                is_active(allocation) ? sum(user_demand.allocated[id.idx, :]) :
+                sum(user_demand.demand[id.idx, :])
+            factor_basin_min * factor_level_min * allocated_total, allocated_total
+        end
         limit_flow!(
             u.user_demand_inflow,
             uprev.user_demand_inflow,
             id,
-            0.0,
-            total_demand,
+            min_flow_rate,
+            max_flow_rate,
             active,
             dt,
         )
     end
 
-    # Evaporation is in [0, ∞) (a stricter upper bound would require also estimating the area)
-    # Infiltration is in [0, infiltration]
+    # Evaporation is in [0, ∞) (stricter bounds would require also estimating the area)
+    # Infiltration is in [f * infiltration, infiltration] where f is a rough estimate of the smallest low storage factor
+    # reduction factor value that was attained over the last timestep
     for (id, infiltration) in zip(basin.node_id, basin.vertical_flux.infiltration)
+        factor_min = min_low_storage_factor(current_storage, basin.storage_prev, id)
         limit_flow!(u.evaporation, uprev.evaporation, id, 0.0, Inf, true, dt)
-        limit_flow!(u.infiltration, uprev.infiltration, id, 0.0, infiltration, true, dt)
+        limit_flow!(
+            u.infiltration,
+            uprev.infiltration,
+            id,
+            factor_min * infiltration,
+            infiltration,
+            true,
+            dt,
+        )
     end
 
     return nothing
diff --git a/core/src/util.jl b/core/src/util.jl
index 77127e99c..d12498c9b 100644
--- a/core/src/util.jl
+++ b/core/src/util.jl
@@ -375,17 +375,11 @@ function reduction_factor(x::T, threshold::Real)::T where {T <: Real}
     end
 end
 
-"If id is a Basin with storage below the threshold, return a reduction factor != 1"
-function low_storage_factor(
-    storage::AbstractVector{T},
+function get_low_storage_factor(
+    current_low_storage_factor::Vector{T},
     id::NodeID,
-    threshold::Real,
-)::T where {T <: Real}
-    if id.type == NodeType.Basin
-        reduction_factor(storage[id.idx], threshold)
-    else
-        one(T)
-    end
+)::T where {T}
+    return id.type == NodeType.Basin ? current_low_storage_factor[id.idx] : one(T)
 end
 
 """
@@ -393,16 +387,15 @@ For resistance nodes, give a reduction factor based on the upstream node
 as defined by the flow direction.
 """
 function low_storage_factor_resistance_node(
-    current_storage,
+    current_low_storage_factor,
     q,
     inflow_id,
     outflow_id,
-    threshold,
 )
     if q > 0
-        low_storage_factor(current_storage, inflow_id, threshold)
+        get_low_storage_factor(current_low_storage_factor, inflow_id)
     else
-        low_storage_factor(current_storage, outflow_id, threshold)
+        get_low_storage_factor(current_low_storage_factor, outflow_id)
     end
 end
 
@@ -648,7 +641,7 @@ function get_variable_ref(
     u = build_state_vector(p)
 
     ref = if node_id.type == NodeType.Basin && variable == "level"
-        PreallocationRef(basin.current_level, node_id.idx)
+        PreallocationRef(basin.current_properties.current_level, node_id.idx)
     elseif variable == "flow_rate" && node_id.type != NodeType.FlowBoundary
         if listen
             if node_id.type ∉ conservative_nodetypes
@@ -861,13 +854,7 @@ end
 
 # Custom overloads
 reduction_factor(x::GradientTracer, threshold::Real) = x
-low_storage_factor_resistance_node(
-    storage,
-    q::GradientTracer,
-    inflow_id,
-    outflow_id,
-    threshold,
-) = q
+low_storage_factor_resistance_node(storage, q::GradientTracer, inflow_id, outflow_id) = q
 relaxed_root(x::GradientTracer, threshold::Real) = x
 get_level_from_storage(basin::Basin, state_idx::Int, storage::GradientTracer) = storage
 stop_declining_negative_storage!(du, u::ComponentVector{<:GradientTracer}) = nothing
@@ -1129,7 +1116,7 @@ end
 Check whether any storages are negative given the state u.
 """
 function isoutofdomain(u, p, t)
-    (; current_storage) = p.basin
+    (; current_storage) = p.basin.current_properties
     current_storage = current_storage[parent(u)]
     formulate_storages!(current_storage, u, u, p, t)
     any(<(0), current_storage)
@@ -1143,3 +1130,44 @@ function get_demand(user_demand, id, priority_idx, t)::Float64
         demand[id.idx, priority_idx]
     end
 end
+
+"""
+Estimate the minimum reduction factor achieved over the last time step by
+estimating the lowest storage achieved over the last time step. To make sure
+it is an underestimate of the minimum, 2LOW_STORAGE_THRESHOLD is subtracted from this lowest storage.
+This is done to not be too strict in clamping the flow in the limiter
+"""
+function min_low_storage_factor(storage_now::Vector{T}, storage_prev, id) where {T}
+    if id.type == NodeType.Basin
+        reduction_factor(
+            min(storage_now[id.idx], storage_prev[id.idx]) - 2LOW_STORAGE_THRESHOLD,
+            LOW_STORAGE_THRESHOLD,
+        )
+    else
+        one(T)
+    end
+end
+
+"""
+Estimate the minimum level reduction factor achieved over the last time step by
+estimating the lowest level achieved over the last time step. To make sure
+it is an underestimate of the minimum, 2USER_DEMAND_MIN_LEVEL_THRESHOLD is subtracted from this lowest level.
+This is done to not be too strict in clamping the flow in the limiter
+"""
+function min_low_user_demand_level_factor(
+    level_now::Vector{T},
+    level_prev,
+    min_level,
+    id_user_demand,
+    id_inflow,
+) where {T}
+    if id_inflow.type == NodeType.Basin
+        reduction_factor(
+            min(level_now[id_inflow.idx], level_prev[id_inflow.idx]) -
+            min_level[id_user_demand.idx] - 2USER_DEMAND_MIN_LEVEL_THRESHOLD,
+            USER_DEMAND_MIN_LEVEL_THRESHOLD,
+        )
+    else
+        one(T)
+    end
+end
diff --git a/core/src/write.jl b/core/src/write.jl
index 425c4beb9..a8ddb5660 100644
--- a/core/src/write.jl
+++ b/core/src/write.jl
@@ -107,7 +107,10 @@ function basin_state_table(
     # ensure the levels are up-to-date
     set_current_basin_properties!(du, u, p, t)
 
-    return (; node_id = Int32.(p.basin.node_id), level = p.basin.current_level[Float64[]])
+    return (;
+        node_id = Int32.(p.basin.node_id),
+        level = p.basin.current_properties.current_level[Float64[]],
+    )
 end
 
 "Create the basin result table from the saved data"
diff --git a/core/test/io_test.jl b/core/test/io_test.jl
index da482360f..36fe71778 100644
--- a/core/test/io_test.jl
+++ b/core/test/io_test.jl
@@ -151,9 +151,10 @@ end
 
     config = Ribasim.Config(toml_path)
     model = Ribasim.Model(config)
-    storage1_begin = copy(model.integrator.p.basin.current_storage[Float64[]])
+    storage1_begin =
+        copy(model.integrator.p.basin.current_properties.current_storage[Float64[]])
     solve!(model)
-    storage1_end = model.integrator.p.basin.current_storage[Float64[]]
+    storage1_end = model.integrator.p.basin.current_properties.current_storage[Float64[]]
     @test storage1_begin != storage1_end
 
     # copy state results to input
@@ -169,6 +170,6 @@ end
     end
 
     model = Ribasim.Model(toml_path)
-    storage2_begin = model.integrator.p.basin.current_storage[Float64[]]
+    storage2_begin = model.integrator.p.basin.current_properties.current_storage[Float64[]]
     @test storage1_end ≈ storage2_begin
 end
diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl
index bfc634ab1..fc6c5e66d 100644
--- a/core/test/run_models_test.jl
+++ b/core/test/run_models_test.jl
@@ -130,7 +130,7 @@ end
     model = Ribasim.run(toml_path)
     @test model isa Ribasim.Model
     (; basin) = model.integrator.p
-    @test basin.current_storage[Float64[]] ≈ [1000]
+    @test basin.current_properties.current_storage[Float64[]] ≈ [1000]
     @test basin.vertical_flux.precipitation == [0.0]
     @test basin.vertical_flux.drainage == [0.0]
     du = get_du(model.integrator)
@@ -153,7 +153,7 @@ end
     du = get_du(integrator)
     (; u, p, t) = integrator
     Ribasim.water_balance!(du, u, p, t)
-    stor = integrator.p.basin.current_storage[parent(du)]
+    stor = integrator.p.basin.current_properties.current_storage[parent(du)]
     prec = p.basin.vertical_flux.precipitation
     evap = du.evaporation
     drng = p.basin.vertical_flux.drainage
@@ -211,7 +211,7 @@ end
 
     @test successful_retcode(model)
     @test length(model.integrator.sol) == 2 # start and end
-    @test model.integrator.p.basin.current_storage[Float64[]] ≈
+    @test model.integrator.p.basin.current_properties.current_storage[Float64[]] ≈
           Float32[828.5386, 801.88289, 492.290, 1318.3053] skip = Sys.isapple() atol = 1.5
 
     @test length(logger.logs) > 10
@@ -248,7 +248,7 @@ end
     du = get_du(model.integrator)
     precipitation = model.integrator.p.basin.vertical_flux.precipitation
     @test length(precipitation) == 4
-    @test model.integrator.p.basin.current_storage[parent(du)] ≈
+    @test model.integrator.p.basin.current_properties.current_storage[parent(du)] ≈
           Float32[721.17656, 695.8066, 416.66188, 1334.4879] atol = 2.0 skip = Sys.isapple()
 end
 
@@ -302,7 +302,7 @@ end
     model = Ribasim.run(toml_path)
     @test model isa Ribasim.Model
     @test successful_retcode(model)
-    @test model.integrator.p.basin.current_storage[Float64[]] ≈
+    @test model.integrator.p.basin.current_properties.current_storage[Float64[]] ≈
           Float32[368.31558, 365.68442] skip = Sys.isapple()
     # the highest level in the dynamic table is updated to 1.2 from the callback
     @test model.integrator.p.tabulated_rating_curve.table[end].t[end] == 1.2
@@ -410,8 +410,9 @@ end
     @test ispath(toml_path)
     model = Ribasim.Model(toml_path)
 
-    (; u, p, t, sol) = model.integrator
-    current_storage = p.basin.current_storage[Float64[]]
+    (; integrator) = model
+    (; u, p, t, sol) = integrator
+    current_storage = p.basin.current_properties.current_storage[Float64[]]
 
     day = 86400.0
 
@@ -505,7 +506,7 @@ end
 
     du = get_du(model.integrator)
     (; p, t) = model.integrator
-    h_actual = p.basin.current_level[parent(du)][1:50]
+    h_actual = p.basin.current_properties.current_level[parent(du)][1:50]
     x = collect(10.0:20.0:990.0)
     h_expected = standard_step_method(x, 5.0, 1.0, 0.04, h_actual[end], 1.0e-6)
 
@@ -592,3 +593,46 @@ end
     @test all(flow .≈ 1.0)
     @test length(flow) == length(tstops) - 1
 end
+
+@testitem "stroboscopic_forcing" begin
+    import BasicModelInterface as BMI
+    using SciMLBase: successful_retcode
+
+    toml_path = normpath(@__DIR__, "../../generated_testmodels/bucket/ribasim.toml")
+    model = BMI.initialize(Ribasim.Model, toml_path)
+
+    drn = BMI.get_value_ptr(model, "basin.drainage")
+    drn_sum = BMI.get_value_ptr(model, "basin.cumulative_drainage")
+    inf = BMI.get_value_ptr(model, "basin.infiltration")
+    inf_sum = BMI.get_value_ptr(model, "basin.cumulative_infiltration")
+
+    nday = 300
+    inf_out = fill(NaN, nday)
+    drn_out = fill(NaN, nday)
+
+    Δt::Float64 = 86400.0
+
+    for day in 0:(nday - 1)
+        if iseven(day)
+            drn .= 25.0 / Δt
+            inf .= 0.0
+        else
+            drn .= 0.0
+            inf .= 25.0 / Δt
+        end
+        BMI.update_until(model, day * Δt)
+
+        inf_out[day + 1] = only(inf_sum)
+        drn_out[day + 1] = only(drn_sum)
+    end
+
+    @test successful_retcode(model)
+
+    Δdrn = diff(drn_out)
+    Δinf = diff(inf_out)
+
+    @test all(Δdrn[1:2:end] .== 0.0)
+    @test all(isapprox.(Δdrn[2:2:end], 25.0; atol = 1e-10))
+    @test all(isapprox.(Δinf[1:2:end], 25.0; atol = 1e-10))
+    @test all(Δinf[2:2:end] .== 0.0)
+end
diff --git a/core/test/utils_test.jl b/core/test/utils_test.jl
index f8a3fd12e..74e082f08 100644
--- a/core/test/utils_test.jl
+++ b/core/test/utils_test.jl
@@ -21,10 +21,6 @@ end
     level_to_area = LinearInterpolation.(area, level)
     storage_to_level = invert_integral.(level_to_area)
     demand = zeros(2)
-    current_level = cache(2)
-    current_area = cache(2)
-    current_level[Float64[]] .= [2.0, 3.0]
-    current_area[Float64[]] .= [2.0, 3.0]
 
     substances = OrderedSet([:test])
     concentration_state = zeros(2, 1)
@@ -33,8 +29,6 @@ end
 
     basin = Ribasim.Basin(;
         node_id = NodeID.(:Basin, [5, 7], [1, 2]),
-        current_level,
-        current_area,
         storage_to_level,
         level_to_area,
         demand,
@@ -46,6 +40,11 @@ end
         concentration_time = StructVector{Ribasim.BasinConcentrationV1}(undef, 0),
     )
 
+    (; current_level, current_area) = basin.current_properties
+
+    current_level[Float64[]] .= [2.0, 3.0]
+    current_area[Float64[]] .= [2.0, 3.0]
+
     @test Ribasim.basin_levels(basin, 2)[1] === 4.0
     @test Ribasim.basin_bottom(basin, NodeID(:Basin, 5, 1))[2] === 0.0
     @test Ribasim.basin_bottom(basin, NodeID(:Basin, 7, 2))[2] === 4.0
@@ -261,28 +260,6 @@ end
     @test reduction_factor(-Inf, 2.0) === 0.0
 end
 
-@testitem "low_storage_factor" begin
-    using Ribasim: NodeID, low_storage_factor, low_storage_factor_resistance_node
-
-    node_id = NodeID(:Basin, 5, 1)
-    @test low_storage_factor([-2.0], node_id, 2.0) === 0.0
-    @test low_storage_factor([0.0f0], node_id, 2.0) === 0.0f0
-    @test low_storage_factor([0.0], node_id, 2.0) === 0.0
-    @test low_storage_factor([1.0f0], node_id, 2.0) === 0.5f0
-    @test low_storage_factor([1.0], node_id, 2.0) === 0.5
-    @test low_storage_factor([3.0f0], node_id, 2.0) === 1.0f0
-    @test low_storage_factor([3.0], node_id, 2.0) === 1.0
-
-    node_id_1 = NodeID(:Basin, 5, 1)
-    node_id_2 = NodeID(:Basin, 6, 2)
-    @test low_storage_factor_resistance_node([3.0, 3.0], 1.0, node_id_1, node_id_2, 2.0) ==
-          1.0
-    @test low_storage_factor_resistance_node([1.0, 3.0], 1.0, node_id_1, node_id_2, 2.0) ==
-          0.5
-    @test low_storage_factor_resistance_node([1.0, 3.0], -1.0, node_id_1, node_id_2, 2.0) ==
-          1.0
-end
-
 @testitem "constraints_from_nodes" begin
     using Ribasim:
         Model,