diff --git a/docs/src/manual/architecture.md b/docs/src/manual/architecture.md
index 7bfd05e41342e286e264f579dc36b951979f825e..eb0720e2b34c9a70854733c0882261822714f621 100644
--- a/docs/src/manual/architecture.md
+++ b/docs/src/manual/architecture.md
@@ -54,7 +54,7 @@ The `CryoGrid` module defines three primary methods that can be used to implemen
 
 1. [`computediagnostic!`](@ref) updates all (non-flux) state variables and/or derived quantities based on the current (prognostic) state.
 2. [`interact!`](@ref) defines interactions between adjacent layers in the stratigraphy, including fluxes over the layer boundary.
-3. [`computefluxes!`](@ref) computes all internal fluxes (and the divergence thereof) within each layer, after boundary fluxes are taken into account by `interact!`.
+3. [`computeprognostic!`](@ref) computes all internal fluxes (and the divergence thereof) within each layer, after boundary fluxes are taken into account by `interact!`.
 
 Layer and/or process specific implementations of each of these methods can generally assume that the previous methods have already been invoked by the caller (it is the responsibility of the calling code to ensure that this is the case). This is, for example, the order in which these methods will be invoked by `tile(du, u, p t)`.
 
@@ -82,7 +82,7 @@ In order to facilitate modularity and ease-of-use, CryoGrid.jl provides an autom
 
 `Prognostic`(@ref) state variables fully define the state of the system at any given time `t`. They form what is typically called the "phase space" or "state space" in the mathematics and engineering literature. In order to be compatible with standard ODE solvers (e.g. like those in `OrdinaryDiffEq`), CryoGrid.jl automatically assembles prognostic state variables into a single array `u` (and its corresponding time derivative `du`) which is returned when initializing a `Tile` with the `initialcondition!` method. Note again that this array should always fully define the state of the system.
 
-`Diagnostic`(@ref) state variables act as caches for intermediate and derived quantities defined by the model. They also may, in some cases, provide a means of coupling between different processes (e.g. the heat and water flux variables `jH` and `jw` might be updated by more than one `Process`). For any model configuration, all diagnostic variables should be fully updated (and thus consistent) with the given prognostic state after invoking `computediagnostic!`, `interact!`, and `computefluxes!`.
+`Diagnostic`(@ref) state variables act as caches for intermediate and derived quantities defined by the model. They also may, in some cases, provide a means of coupling between different processes (e.g. the heat and water flux variables `jH` and `jw` might be updated by more than one `Process`). For any model configuration, all diagnostic variables should be fully updated (and thus consistent) with the given prognostic state after invoking `computediagnostic!`, `interact!`, and `computeprognostic!`.
 
 When a `Tile` is constructed, all variables defined by each layer in the `Stratigraphy` are collected and then intiailized in [`StateVars`](@ref) according to the given `DiscretizationStrategy`.
 
diff --git a/docs/src/manual/overview.md b/docs/src/manual/overview.md
index bb5cd1ee25704ccceac96d5ed211f58273473b88..e58d8030fc029780948a000c2996cfdc20e9425e 100644
--- a/docs/src/manual/overview.md
+++ b/docs/src/manual/overview.md
@@ -69,12 +69,12 @@ When the `HeatBalance` process is assigned to a `Soil` layer, `Tile` will invoke
 
 Each variable definition consists of a name (a Julia `Symbol`), a type, and a shape. For variables discretized on the grid, the shape is specified by `OnGrid`, which will generate an array of the appropriate size when the model is compiled. The arguments `Cells` and `Edges` specify whether the variable should be defined on the grid cells or edges respecitvely.
 
-The real work finally happens in [`computediagnostic!`](@ref) and [`computefluxes!`](@ref), the latter of which should be used to compute the time derivatives (here `dH`). [`interact!`](@ref) defines the behavior at the boundaries and should be used to compute the derivatives (and any other necessary values) at the interface between layers.
+The real work finally happens in [`computediagnostic!`](@ref) and [`computeprognostic!`](@ref), the latter of which should be used to compute the time derivatives (here `dH`). [`interact!`](@ref) defines the behavior at the boundaries and should be used to compute the derivatives (and any other necessary values) at the interface between layers.
 
-We can take as an example the implementation of `computefluxes!` for enthalpy-based heat conduction (note that `jH` is a diagnostic variable representing the energy flux over each cell edge):
+We can take as an example the implementation of `computeprognostic!` for enthalpy-based heat conduction (note that `jH` is a diagnostic variable representing the energy flux over each cell edge):
 
 ```julia
-function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:EnthalpyBased}, state)
+function CryoGrid.computeprognostic!(::SubSurface, ::HeatBalance{<:EnthalpyBased}, state)
     Δk = Δ(state.grid) # cell sizes
     ΔT = Δ(cells(state.grid)) # midpoint distances
     # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
diff --git a/src/CryoGrid.jl b/src/CryoGrid.jl
index 97c9a9944a1dbb63027ea5824ee4907d09d247d9..c3dfbf7814a95f2e865d111a4c5c0fc24719c1cb 100755
--- a/src/CryoGrid.jl
+++ b/src/CryoGrid.jl
@@ -54,7 +54,7 @@ include("variables.jl")
 export BCKind
 include("traits.jl")
 
-export initialcondition!, computediagnostic!, interact!, interactmaybe!, computefluxes!, resetfluxes!, diagnosticstep!
+export initialcondition!, computediagnostic!, interact!, interactmaybe!, computeprognostic!, resetfluxes!, diagnosticstep!
 export variables, processes, initializers, timestep, isactive, caninteract
 export boundaryflux, boundaryvalue, criterion, criterion!, trigger!
 include("methods.jl")
diff --git a/src/Physics/Heat/heat_conduction.jl b/src/Physics/Heat/heat_conduction.jl
index 6a237d616a506920d01035c50434300c7a1b03a6..785230f74da2649b4968391a4659bec16745eee8 100644
--- a/src/Physics/Heat/heat_conduction.jl
+++ b/src/Physics/Heat/heat_conduction.jl
@@ -104,7 +104,7 @@ function CryoGrid.interact!(sub1::SubSurface, ::HeatBalance, sub2::SubSurface, :
     return nothing
 end
 
-function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:EnthalpyBased}, state)
+function CryoGrid.computeprognostic!(::SubSurface, ::HeatBalance{<:EnthalpyBased}, state)
     Δk = Δ(state.grid) # cell sizes
     ΔT = Δ(cells(state.grid)) # midpoint distances
     # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set;
@@ -115,7 +115,7 @@ function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:EnthalpyBased}, s
     divergence!(state.dH, state.jH, Δk)
     return nothing
 end
-function CryoGrid.computefluxes!(sub::SubSurface, ::HeatBalance{<:TemperatureBased}, state)
+function CryoGrid.computeprognostic!(sub::SubSurface, ::HeatBalance{<:TemperatureBased}, state)
     Δk = Δ(state.grid) # cell sizes
     ΔT = Δ(cells(state.grid)) # midpoint distances
     # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
diff --git a/src/Physics/Heat/heat_implicit.jl b/src/Physics/Heat/heat_implicit.jl
index dded8b6fbc456509b4459b8d3c42386fbded221e..254527d8d0e65c4c9aa831ed941375896a9e3d60 100644
--- a/src/Physics/Heat/heat_implicit.jl
+++ b/src/Physics/Heat/heat_implicit.jl
@@ -82,8 +82,8 @@ function CryoGrid.interact!(sub1::SubSurface, ::HeatBalanceImplicit, sub2::SubSu
     return nothing
 end
 
-# do nothing in computefluxes!
-CryoGrid.computefluxes!(::SubSurface, ::HeatBalanceImplicit, state) = nothing
+# do nothing in computeprognostic!
+CryoGrid.computeprognostic!(::SubSurface, ::HeatBalanceImplicit, state) = nothing
 
 function CryoGrid.resetfluxes!(sub::SubSurface, heat::HeatBalanceImplicit, state)
     @inbounds for i in 1:length(state.H)
diff --git a/src/Physics/Heat/heat_water.jl b/src/Physics/Heat/heat_water.jl
index 5c513182411145961fb8d90c2f93e4fcab1efe60..8ea5878298ba8d9961c3055864e05acdef3b11d2 100644
--- a/src/Physics/Heat/heat_water.jl
+++ b/src/Physics/Heat/heat_water.jl
@@ -89,13 +89,13 @@ function CryoGrid.interact!(
 end
 
 # Flux calculation
-function CryoGrid.computefluxes!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
+function CryoGrid.computeprognostic!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
     water, heat = ps
-    CryoGrid.computefluxes!(sub, water, state)
+    CryoGrid.computeprognostic!(sub, water, state)
     if heat.advection
         water_energy_advection!(sub, ps, state)
     end
-    CryoGrid.computefluxes!(sub, heat, state)
+    CryoGrid.computeprognostic!(sub, heat, state)
 end
 
 function heatcapacitywater(sub::SubSurface, state)
diff --git a/src/Physics/Hydrology/water_balance.jl b/src/Physics/Hydrology/water_balance.jl
index b87ff3ef4ba37a0c7264b86735f6649f214a58a6..0aea1b9ec3c44bfe39a4f73faaab3cf55f5d4eb0 100644
--- a/src/Physics/Hydrology/water_balance.jl
+++ b/src/Physics/Hydrology/water_balance.jl
@@ -219,7 +219,7 @@ function CryoGrid.computediagnostic!(sub::SubSurface, water::WaterBalance, state
     evapotranspiration!(sub, water, state)
 end
 
-function CryoGrid.computefluxes!(sub::SubSurface, water::WaterBalance, state)
+function CryoGrid.computeprognostic!(sub::SubSurface, water::WaterBalance, state)
     wateradvection!(sub, water, state)
     waterdiffusion!(sub, water, state)
     balancefluxes!(sub, water, state)
@@ -255,7 +255,7 @@ function CryoGrid.initialcondition!(sub::SubSurface, water::WaterBalance{NoFlow}
     end
 end
 
-CryoGrid.computefluxes!(::SubSurface, ::WaterBalance{NoFlow}, state) = nothing
+CryoGrid.computeprognostic!(::SubSurface, ::WaterBalance{NoFlow}, state) = nothing
 
 CryoGrid.resetfluxes!(::SubSurface, ::WaterBalance{NoFlow}, state) = nothing
 
diff --git a/src/Physics/Lakes/lake_simple.jl b/src/Physics/Lakes/lake_simple.jl
index 643e887edf75f1ce52350092f8f3b3d7ee172a04..7a1ae48dedf522bf98d473bcbbcad49cc2a09e80 100644
--- a/src/Physics/Lakes/lake_simple.jl
+++ b/src/Physics/Lakes/lake_simple.jl
@@ -110,7 +110,7 @@ function CryoGrid.computediagnostic!(sub::Lake, heat::HeatBalance{<:Heat.Diffusi
     return nothing
 end
 
-function CryoGrid.computefluxes!(::Lake, ::HeatBalance{<:Heat.Diffusion1D{:H}}, state)
+function CryoGrid.computeprognostic!(::Lake, ::HeatBalance{<:Heat.Diffusion1D{:H}}, state)
     Δk = Δ(state.grids.k) # cell sizes
     ΔT = Δ(state.grids.T) # midpoint distances
     # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
diff --git a/src/Physics/Salt/salt_diffusion.jl b/src/Physics/Salt/salt_diffusion.jl
index f6175a76e6e2ca37356619d5e28e3d1f6f4f97c8..1401acc4c1109365afef7eba8af88967e5f19f36 100644
--- a/src/Physics/Salt/salt_diffusion.jl
+++ b/src/Physics/Salt/salt_diffusion.jl
@@ -121,7 +121,7 @@ function CryoGrid.timestep(::SalineGround, salt::SaltMassBalance{T,<:CryoGrid.CF
     return dtmax
 end
 
-function CryoGrid.computefluxes!(
+function CryoGrid.computeprognostic!(
     ::SalineGround,
     ps::CoupledHeatSalt{THeat},
     state
diff --git a/src/Physics/Snow/snow_bulk.jl b/src/Physics/Snow/snow_bulk.jl
index 230a7efa4094c219f0d6cfd475454278c202bd65..35a7cbdff37d45f87720db5f60ed9362482e2291 100644
--- a/src/Physics/Snow/snow_bulk.jl
+++ b/src/Physics/Snow/snow_bulk.jl
@@ -167,16 +167,16 @@ function CryoGrid.trigger!(
 end
 
 # mass balance fluxes are handled in interact! for bulk snowpack
-CryoGrid.computefluxes!(snow::BulkSnowpack, mass::SnowMassBalance, state) = nothing
+CryoGrid.computeprognostic!(snow::BulkSnowpack, mass::SnowMassBalance, state) = nothing
 
-# computefluxes! for free water, enthalpy based HeatBalance on bulk snow layer
-function CryoGrid.computefluxes!(
+# computeprognostic! for free water, enthalpy based HeatBalance on bulk snow layer
+function CryoGrid.computeprognostic!(
     snow::BulkSnowpack,
     ps::CoupledSnowWaterHeat{TM,TW,TH},
     state
 ) where {TM,TW,TH<:HeatBalance{<:EnthalpyBased}}
     mass, water, heat = ps
-    computefluxes!(snow, mass, state)
+    computeprognostic!(snow, mass, state)
     dsn = getscalar(state.dsn)
     if dsn < snow.para.thresh
         # set fluxes to zero if there is no snow
@@ -185,8 +185,8 @@ function CryoGrid.computefluxes!(
             state.dsat .= 0.0
         end
     else
-        # otherwise call computefluxes! for coupled water/heat
-        computefluxes!(snow, Coupled(water, heat), state)
+        # otherwise call computeprognostic! for coupled water/heat
+        computeprognostic!(snow, Coupled(water, heat), state)
     end
     return nothing
 end
diff --git a/src/Physics/Snow/snow_lite.jl b/src/Physics/Snow/snow_lite.jl
index 597194d3b6bd3de1fc8c098c2a955cc075cc61e1..1dbbeb700772d977bd16e72c5615a427e32c468b 100644
--- a/src/Physics/Snow/snow_lite.jl
+++ b/src/Physics/Snow/snow_lite.jl
@@ -40,7 +40,7 @@ CryoGrid.variables(::LiteSnowpack, ::SnowMassBalance) = (
 
 CryoGrid.computediagnostic!(::LiteSnowpack, ::SnowMassBalance, state) = nothing
 
-CryoGrid.computefluxes!(::LiteSnowpack, ::SnowMassBalance, state) = nothing
+CryoGrid.computeprognostic!(::LiteSnowpack, ::SnowMassBalance, state) = nothing
 
 function CryoGrid.interact!(
     top::Top,
diff --git a/src/Physics/Snow/snowpack.jl b/src/Physics/Snow/snowpack.jl
index 69d6705044e2c79dd5e91a98c0af3951c40637fc..4644619cf28796433e1be0f01465b73f15c1e5c0 100644
--- a/src/Physics/Snow/snowpack.jl
+++ b/src/Physics/Snow/snowpack.jl
@@ -120,14 +120,14 @@ function CryoGrid.computediagnostic!(
     updategrid!(state.grid, z0, state.dsn)
 end
 
-function CryoGrid.computefluxes!(
+function CryoGrid.computeprognostic!(
     snow::Snowpack,
     state,
 )
     mass, water, heat = processes(snow)
-    computefluxes!(snow, mass, state)
-    computefluxes!(snow, water, state)
-    computefluxes!(snow, heat, state)
+    computeprognostic!(snow, mass, state)
+    computeprognostic!(snow, water, state)
+    computeprognostic!(snow, heat, state)
 end
 
 # Special overrides for heat timestep control on snow layer
diff --git a/src/Physics/Sources/Sources.jl b/src/Physics/Sources/Sources.jl
index 20baf9b080372adf8020d96ce560319ad963dff1..53b58d4dd1917fcf89c087c8b67ad2a7863f30dc 100644
--- a/src/Physics/Sources/Sources.jl
+++ b/src/Physics/Sources/Sources.jl
@@ -1,7 +1,7 @@
 module Sources
 
 import CryoGrid: SubSurfaceProcess, SubSurface
-import CryoGrid: computediagnostic!, initialcondition!, interact!, computefluxes!, variables
+import CryoGrid: computediagnostic!, initialcondition!, interact!, computeprognostic!, variables
 
 using ..Heat
 using CryoGrid.Numerics
@@ -52,7 +52,7 @@ ConstructionBase.constructorof(::Type{<:Source{P}}) where {P} = (term, aux) -> S
 
 (p::Periodic)(t) = p.amp*sin(2Ï€*p.freq*t - p.shift) + p.level
 # HeatBalance sources
-computefluxes!(::SubSurface, s::Source{<:HeatBalance,<:Constant}, state) = @inbounds @. state.dH += s.term.Sâ‚€
-computefluxes!(::SubSurface, src::Source{<:HeatBalance,<:Periodic}, state) = @inbounds @. state.dH += src.term(state.t)
+computeprognostic!(::SubSurface, s::Source{<:HeatBalance,<:Constant}, state) = @inbounds @. state.dH += s.term.Sâ‚€
+computeprognostic!(::SubSurface, src::Source{<:HeatBalance,<:Periodic}, state) = @inbounds @. state.dH += src.term(state.t)
 
 end
\ No newline at end of file
diff --git a/src/Physics/Surface/SWB/coupled_sweb.jl b/src/Physics/Surface/SWB/coupled_sweb.jl
index 19753844faa3ab37972afad61e3934e80e2c0bda..53a5de7929f10f3ca873209d6af5eef85fccadd8 100644
--- a/src/Physics/Surface/SWB/coupled_sweb.jl
+++ b/src/Physics/Surface/SWB/coupled_sweb.jl
@@ -13,11 +13,11 @@ CryoGrid.variables(top::Top, bc::SurfaceWaterEnergyBalance) = (
     Prognostic(:ET, Scalar, u"m^3"),
 )
 
-function CryoGrid.computefluxes!(top::Top, bc::SurfaceWaterEnergyBalance, state)
+function CryoGrid.computeprognostic!(top::Top, bc::SurfaceWaterEnergyBalance, state)
     swb, seb = bc
-    computefluxes!(top, swb, state)
+    computeprognostic!(top, swb, state)
     ET!(top, swb, state)
-    computefluxes!(top, seb, state)
+    computeprognostic!(top, seb, state)
 end
 
 function CryoGrid.interact!(top::Top, bc::SurfaceWaterEnergyBalance, sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), stop, ssub)
diff --git a/src/Physics/Surface/SWB/swb.jl b/src/Physics/Surface/SWB/swb.jl
index c395e12960807a03152d66fed8fd824e21a00391..3003d3df44e0e77701f34c0d2454533d49d77c17 100644
--- a/src/Physics/Surface/SWB/swb.jl
+++ b/src/Physics/Surface/SWB/swb.jl
@@ -71,7 +71,7 @@ function CryoGrid.computediagnostic!(::Top, swb::SurfaceWaterBalance, stop)
     @setscalar stop.jw_rain = swb.rainfall(stop.t)
 end
 
-function CryoGrid.computefluxes!(top::Top, swb::SurfaceWaterBalance, state)
+function CryoGrid.computeprognostic!(top::Top, swb::SurfaceWaterBalance, state)
     runoff!(top, swb, state)
 end
 
diff --git a/src/Tiles/stratigraphy.jl b/src/Tiles/stratigraphy.jl
index a54157c958ba206314aa5741ec19335cd5ef6670..e7fab3649dbdfbb69288c3cd727bdae78a642c72 100644
--- a/src/Tiles/stratigraphy.jl
+++ b/src/Tiles/stratigraphy.jl
@@ -219,9 +219,9 @@ that `getproperty` would return the appropriate state object for the i'th layer
     return expr
 end
 
-function CryoGrid.computefluxes!(strat::Stratigraphy, state)
+function CryoGrid.computeprognostic!(strat::Stratigraphy, state)
     fastiterate(namedlayers(strat)) do named_layer
-        CryoGrid.computefluxes!(named_layer.val, getproperty(state, nameof(named_layer)))
+        CryoGrid.computeprognostic!(named_layer.val, getproperty(state, nameof(named_layer)))
     end
 end
 
diff --git a/src/Tiles/tile.jl b/src/Tiles/tile.jl
index 2b8975591a83e5e10c8747b17bb5a777ac4c3ae4..ed4d2c295b07729492f45188d08f55e24c2bd703 100644
--- a/src/Tiles/tile.jl
+++ b/src/Tiles/tile.jl
@@ -108,7 +108,7 @@ function Tiles.Tile(integrator::SciMLBase.DEIntegrator)
 end
 
 """
-    computefluxes!(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents}
+    computeprognostic!(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents}
 
 Time derivative step function (i.e. du/dt) for any arbitrary `Tile`. Specialized code is generated and compiled
 on the fly via the @generated macro to ensure type stability. The generated code updates each layer in the stratigraphy
@@ -117,10 +117,10 @@ in sequence, i.e for each layer 1 <= i <= N:
 ```julia
 computediagnostic!(layer[i], ...)
 interact!(layer[i], ..., layer[i+1], ...)
-computefluxes!(layer[i], ...)
+computeprognostic!(layer[i], ...)
 ```
 """
-function computefluxes!(
+function computeprognostic!(
     _tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true},
     _du::AbstractVector,
     _u::AbstractVector,
@@ -139,8 +139,8 @@ function computefluxes!(
     checkstate!(tile, state, u, du, :computediagnostic!)
     CryoGrid.interact!(strat, state)
     checkstate!(tile, state, u, du, :interact!)
-    CryoGrid.computefluxes!(strat, state)
-    checkstate!(tile, state, u, du, :computefluxes!)
+    CryoGrid.computeprognostic!(strat, state)
+    checkstate!(tile, state, u, du, :computeprognostic!)
     return nothing
 end
 
@@ -182,7 +182,7 @@ function CryoGrid.initialcondition!(tile::Tile, tspan::NTuple{2,Float64}, p::Abs
     CryoGrid.initialcondition!(tile.grid, state)
     CryoGrid.initialcondition!(strat, state, tile.inits...)
     # evaluate initial time derivative
-    computefluxes!(tile, du, u, p, t0)
+    computeprognostic!(tile, du, u, p, t0)
     return u, du
 end
 
diff --git a/src/Tiles/tile_base.jl b/src/Tiles/tile_base.jl
index 8340d52affaaa178b0ac422ec65ccd3a925a085c..be02dd9149de9c3c98d1a0b9ee9aa2049bf3bce2 100644
--- a/src/Tiles/tile_base.jl
+++ b/src/Tiles/tile_base.jl
@@ -16,9 +16,9 @@ abstract type AbstractTile{iip} end
     (tile::AbstractTile{true})(du,u,p,t,dt=1.0)
     (tile::AbstractTile{false})(u,p,t,dt=1.0)
 
-Invokes `computefluxes!` on `tile` to compute the time derivative du/dt.
+Invokes `computeprognostic!` on `tile` to compute the time derivative du/dt.
 """
-(tile::AbstractTile{true})(du, u, p, t, dt=1.0) = computefluxes!(tile,du,u,p,t,dt)
+(tile::AbstractTile{true})(du, u, p, t, dt=1.0) = computeprognostic!(tile,du,u,p,t,dt)
 (tile::AbstractTile{false})(u, p, t, dt=1.0) = computefluxes(tile,u,p,t,dt)
 
 """
@@ -29,11 +29,11 @@ Returns true if `tile` uses in-place mode, false if out-of-place.
 isinplace(::AbstractTile{iip}) where {iip} = iip
 
 """
-    computefluxes!(::T,du,u,p,t) where {T<:AbstractTile}
+    computeprognostic!(::T,du,u,p,t) where {T<:AbstractTile}
 
 In-place update function for tile `T`. Computes du/dt and stores the result in `du`.
 """
-computefluxes!(::T, du, u, p, t, dt=1.0) where {T<:AbstractTile} = error("no implementation of in-place computefluxes! for $T")
+computeprognostic!(::T, du, u, p, t, dt=1.0) where {T<:AbstractTile} = error("no implementation of in-place computeprognostic! for $T")
 
 """
     computefluxes(::T,u,p,t) where {T<:AbstractTile}
diff --git a/src/coupling.jl b/src/coupling.jl
index 1567f3b0cc0cd7aebb683717b5e49b8980303f1c..785975a68005d150dc9593457ada3a08f6d569f7 100644
--- a/src/coupling.jl
+++ b/src/coupling.jl
@@ -17,7 +17,7 @@ CryoGrid.initialcondition!(layer::Layer, ps::CoupledProcesses, state) = _invoke_
 
 CryoGrid.computediagnostic!(layer::Layer, ps::CoupledProcesses, state) = _invoke_sequential(computediagnostic!, layer, ps, state)
 
-CryoGrid.computefluxes!(layer::Layer, ps::CoupledProcesses, state) = _invoke_sequential(computefluxes!, layer, ps, state)
+CryoGrid.computeprognostic!(layer::Layer, ps::CoupledProcesses, state) = _invoke_sequential(computeprognostic!, layer, ps, state)
 
 CryoGrid.resetfluxes!(layer::Layer, ps::CoupledProcesses, state) = _invoke_sequential(resetfluxes!, layer, ps, state)
 
diff --git a/src/methods.jl b/src/methods.jl
index a78dce27ee8e3947d0d6f4d8151265bb429caaed..b9511218b2698fe8532055d60e2de9327cf7252c 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -63,15 +63,15 @@ interact!(layer1::Layer, layer2::Layer, state1, state2) = interact!(layer1, proc
 interact!(::Layer, ::Process, ::Layer, ::Process, state1, state2) = nothing
 
 """
-    computefluxes!(l::Layer, p::Process, state)
+    computeprognostic!(l::Layer, p::Process, state)
 
-Calculates all internal fluxes for a given layer. Note that an instance of `computefluxes!` must be provided
+Calculates all internal fluxes for a given layer. Note that an instance of `computeprognostic!` must be provided
 for all non-boundary (subsurface) processes/layers.
 """
-computefluxes!(layer::Layer, state) = computefluxes!(layer, processes(layer), state)
-computefluxes!(layer::Layer, proc::Process, state) = error("computefluxes! defined for $(typeof(layer)) with $(typeof(proc))")
-computefluxes!(::Top, ::BoundaryProcess, state) = nothing
-computefluxes!(::Bottom, ::BoundaryProcess, state) = nothing
+computeprognostic!(layer::Layer, state) = computeprognostic!(layer, processes(layer), state)
+computeprognostic!(layer::Layer, proc::Process, state) = error("computeprognostic! defined for $(typeof(layer)) with $(typeof(proc))")
+computeprognostic!(::Top, ::BoundaryProcess, state) = nothing
+computeprognostic!(::Bottom, ::BoundaryProcess, state) = nothing
 
 """
     caninteract(layer1::Layer, layer2::Layer, state1, state2)
@@ -111,7 +111,7 @@ resetfluxes!(layer::Layer, proc::Process, state) = nothing
     isactive(::Layer, state)
 
 Returns a boolean whether or not this layer is currently active in the stratigraphy and should interact with other layers.
-Note that `computediagnostic!` and `computefluxes!` are always invoked regardless of the current state of `isactive`.
+Note that `computediagnostic!` and `computeprognostic!` are always invoked regardless of the current state of `isactive`.
 The default implementation of `isactive` always returns `true`.
 """
 isactive(::Layer, state) = true
diff --git a/test/Physics/Heat/heat_conduction_tests.jl b/test/Physics/Heat/heat_conduction_tests.jl
index 9a2df9563ad30471604a55db7c5e7eea160ee980..692cf824dfe689a8d26e7b50402962f120be151c 100644
--- a/test/Physics/Heat/heat_conduction_tests.jl
+++ b/test/Physics/Heat/heat_conduction_tests.jl
@@ -117,7 +117,7 @@ end
 		state = (T=T,dH=dH,jH=jH,k=k,grid=x,grids=(T=xc,k=x),t=t)
 		interact!(Top(bc), bc, sub, heat, state, state)
 		interact!(sub, heat, Bottom(bc), bc, state, state)
-		computefluxes!(sub, heat, state)
+		computeprognostic!(sub, heat, state)
 		# strip units from dH before returning it to the solver;
 		# note that we do not need to divide by diffusivity since we assume it to be unity
 		return ustrip.(dH)
diff --git a/test/Physics/Hydrology/water_balance_tests.jl b/test/Physics/Hydrology/water_balance_tests.jl
index 2c10d4f39ea397f6347d2eb8baee7a10a2eec591..ac819ccd68c33f4b23fc9c577d62644548124063 100644
--- a/test/Physics/Hydrology/water_balance_tests.jl
+++ b/test/Physics/Hydrology/water_balance_tests.jl
@@ -32,14 +32,14 @@ using Test
         @test allfinite(state.kw)
         @test allfinite(state.θw)
     end
-    @testset "computefluxes!" begin
+    @testset "computeprognostic!" begin
         sub = TestGroundLayer(WaterBalance(BucketScheme()))
         state = Diagnostics.build_dummy_state(testgrid, sub, with_units=true)
         # initialize fully saturated
         state.sat .= 1.0
         procs = CryoGrid.processes(sub)
         CryoGrid.computediagnostic!(sub, procs, state)
-        CryoGrid.computefluxes!(sub, procs, state)
+        CryoGrid.computeprognostic!(sub, procs, state)
         @test allfinite(state.kw)
         @test allfinite(state.θw)
         @test all(iszero.(state.jw))
@@ -47,7 +47,7 @@ using Test
         state.sat .= 0.5
         procs = CryoGrid.processes(sub)
         CryoGrid.computediagnostic!(sub, procs, state)
-        CryoGrid.computefluxes!(sub, procs, state)
+        CryoGrid.computeprognostic!(sub, procs, state)
         @test allfinite(state.kw)
         @test allfinite(state.θw)
         @test all(state.jw[2:end-1] .> zero(eltype(state.jw)))
diff --git a/test/Physics/Salt/runtests.jl b/test/Physics/Salt/runtests.jl
index 8363c46cef3cb6044127ca31d8ce226b87910287..d9dc484526d0bfc2b434dbe8a266f81a52cbd855 100644
--- a/test/Physics/Salt/runtests.jl
+++ b/test/Physics/Salt/runtests.jl
@@ -63,7 +63,7 @@ using Test
         @test state1.jc[end] > zero(eltype(state1.jc))
         @test state2.jc[1] == state1.jc[end]
     end
-    @testset "computefluxes!" begin
+    @testset "computeprognostic!" begin
         soil = pstrip(SalineGround())
         state = Diagnostics.build_dummy_state(testgrid, soil, with_units=false)
         # initialize variables
@@ -72,7 +72,7 @@ using Test
         state.por .= porosity(soil)
         procs = CryoGrid.processes(soil)
         CryoGrid.computediagnostic!(soil, procs, state)
-        CryoGrid.computefluxes!(soil, procs, state)
+        CryoGrid.computeprognostic!(soil, procs, state)
         # check that all divergences are nonzero for both temperature and salt
         @test all(.!iszero.(state.dc))
         @test all(.!iszero.(state.dT))
diff --git a/test/Physics/Sources/runtests.jl b/test/Physics/Sources/runtests.jl
index 1338b5c218afc3c60cb6bc71ca1922d3267f2b50..0b72ea448db5b4bf6ced1f54469ebc8fec5284ff 100644
--- a/test/Physics/Sources/runtests.jl
+++ b/test/Physics/Sources/runtests.jl
@@ -6,10 +6,10 @@ using Test
         heatsource = Source(HeatBalance, Sources.Constant(1.0u"W/m^3"))
         layer = TestGroundLayer(heatsource)
         state = (dH=zeros(100)u"W/m^3",)
-        computefluxes!(layer, heatsource, state)
+        computeprognostic!(layer, heatsource, state)
         @test all(state.dH .≈ 1.0u"W/m^3")
         state = (dH=ones(100)u"W/m^3",)
-        computefluxes!(layer, heatsource, state)
+        computeprognostic!(layer, heatsource, state)
         @test all(state.dH .≈ 2.0u"W/m^3")
     end
     @testset "Periodic" begin
@@ -17,10 +17,10 @@ using Test
         heatsource = Source(HeatBalance, Sources.Periodic(; level, amp, freq, shift))
         layer = TestGroundLayer(heatsource)
         state = (dH=zeros(100)u"W/m^3", t=1.0u"s")
-        computefluxes!(layer, heatsource, state)
+        computeprognostic!(layer, heatsource, state)
         @test all(state.dH .≈ level + amp*sin(2π*freq*1.0u"s" - shift))
         state = (dH=ones(100)u"W/m^3", t=1.5u"s")
-        computefluxes!(layer, heatsource, state)
+        computeprognostic!(layer, heatsource, state)
         @test all(state.dH .≈ 1.0u"W/m^3" + level + amp*sin(2π*freq*1.5u"s" - shift))
     end
 end
diff --git a/test/Physics/Surface/swb_tests.jl b/test/Physics/Surface/swb_tests.jl
index c6aabd6b7def090269656a712af194f915830494..6a52498477fa5069a69fb840d0cf33fca8ef900d 100644
--- a/test/Physics/Surface/swb_tests.jl
+++ b/test/Physics/Surface/swb_tests.jl
@@ -17,7 +17,7 @@
         CryoGrid.computediagnostic!(top, swb, stop)
         CryoGrid.computediagnostic!(sub, water, ssub)
         CryoGrid.interact!(top, sub, stop, ssub)
-        CryoGrid.computefluxes!(top, stop)
+        CryoGrid.computeprognostic!(top, stop)
         @test iszero(ssub.jw[1])
         @test stop.drunoff[1] > zero(eltype(stop.drunoff))
         ssub.sat .= 0.5
@@ -26,7 +26,7 @@
         CryoGrid.computediagnostic!(top, swb, stop)
         CryoGrid.computediagnostic!(sub, water, ssub)
         CryoGrid.interact!(top, sub, stop, ssub)
-        CryoGrid.computefluxes!(top, stop)
+        CryoGrid.computeprognostic!(top, stop)
         @test ssub.jw[1] > zero(eltype(ssub.jw))
         @test iszero(stop.drunoff[1])
     end