diff --git a/Project.toml b/Project.toml
index 2e574ea..f61583c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,9 +15,10 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [targets]
-test = ["Test"]
+test = ["Test", "Statistics"]
 
 [compat]
 CairoMakie = "0.14.0"
diff --git a/test/regression/cyclic_regression.jl b/test/regression/cyclic_regression.jl
new file mode 100644
index 0000000..71a716b
--- /dev/null
+++ b/test/regression/cyclic_regression.jl
@@ -0,0 +1,90 @@
+using Test
+using Mocca
+using Jutul
+using Statistics: mean
+
+function run_cyclic_simulation(; ncells = 200, num_cycles = 3, max_dt = 1.0)
+    constants = Mocca.HaghpanahConstants{Float64}()
+    system = Mocca.TwoComponentAdsorptionSystem(constants)
+    model = Mocca.setup_process_model(system; ncells = ncells)
+
+    bar = si_unit(:bar)
+    state0 = Mocca.setup_process_state(model;
+        Pressure = 1 * bar,
+        Temperature = 298.15,
+        WallTemperature = constants.T_a,
+        y = [1e-10, 1.0 - 1e-10]
+    )
+
+    parameters = Mocca.setup_process_parameters(model)
+
+    stage_times = [15.0, 15.0, 30.0, 40.0]
+    stage_names = ["pressurisation", "adsorption", "blowdown", "evacuation"]
+    sim_forces, timesteps = Mocca.setup_forces(model, stage_times, stage_names;
+        num_cycles = num_cycles, max_dt = max_dt)
+
+    case = Mocca.MoccaCase(model, timesteps, sim_forces;
+        state0 = state0, parameters = parameters)
+
+    states, timesteps_out = Mocca.simulate_process(case;
+        output_substates = true,
+        info_level = -1
+    )
+
+    return states, model, timesteps_out, num_cycles, stage_times
+end
+
+function cyclic_metrics(states, model, timesteps_out, num_cycles, stage_times)
+    nc = number_of_cells(model.domain)
+    final = states[end]
+
+    # Find start of last cycle by tracking cumulative time
+    cycle_duration = sum(stage_times)
+    last_cycle_start_time = (num_cycles - 1) * cycle_duration
+
+    # Find index where last cycle starts
+    cumtime = 0.0
+    last_cycle_start_idx = 1
+    for i in 1:length(timesteps_out)
+        if cumtime >= last_cycle_start_time
+            last_cycle_start_idx = i
+            break
+        end
+        cumtime += timesteps_out[i]
+    end
+
+    last_cycle_states = states[last_cycle_start_idx:end]
+
+    return (
+        y_CO2_final = final[:y][1, nc],
+        y_N2_final = final[:y][2, nc],
+        T_final = final[:Temperature][nc],
+        P_final = final[:Pressure][nc],
+        total_CO2_adsorbed_final = sum(final[:AdsorbedConcentration][1, :]),
+        avg_T_last_cycle = mean([maximum(s[:Temperature]) for s in last_cycle_states]),
+        avg_P_last_cycle = mean([mean(s[:Pressure]) for s in last_cycle_states]),
+    )
+end
+
+@testset "Cyclic VSA Regression" begin
+    states, model, timesteps_out, num_cycles, stage_times = run_cyclic_simulation()
+    m = cyclic_metrics(states, model, timesteps_out, num_cycles, stage_times)
+
+    ref = (
+        y_CO2_final = 7.100524006314865e-8,
+        y_N2_final = 0.9999999289947609,
+        T_final = 295.76741982625623,
+        P_final = 9998.130045421889,
+        total_CO2_adsorbed_final = 52325.83782488843,
+        avg_T_last_cycle = 348.8714798098076,
+        avg_P_last_cycle = 40071.222347007766,
+    )
+
+    @test isapprox(m.y_CO2_final, ref.y_CO2_final, rtol=1e-3, atol=1e-10)
+    @test isapprox(m.y_N2_final, ref.y_N2_final, rtol=1e-3)
+    @test isapprox(m.T_final, ref.T_final, atol=0.5)
+    @test isapprox(m.P_final, ref.P_final, rtol=1e-3)
+    @test isapprox(m.total_CO2_adsorbed_final, ref.total_CO2_adsorbed_final, rtol=1e-3)
+    @test isapprox(m.avg_T_last_cycle, ref.avg_T_last_cycle, atol=0.5)
+    @test isapprox(m.avg_P_last_cycle, ref.avg_P_last_cycle, rtol=1e-3)
+end
diff --git a/test/regression/dcb_regression.jl b/test/regression/dcb_regression.jl
new file mode 100644
index 0000000..d3e6079
--- /dev/null
+++ b/test/regression/dcb_regression.jl
@@ -0,0 +1,86 @@
+using Test
+using Mocca
+using Jutul
+
+function run_dcb_simulation(; ncells = 200, t_ads = 5000.0, maxdt = 5000.0)
+    constants = Mocca.HaghpanahConstants{Float64}(h_in = 0.0, h_out = 0.0)
+    system = Mocca.TwoComponentAdsorptionSystem(constants)
+    model = Mocca.setup_process_model(system; ncells = ncells)
+
+    bar = si_unit(:bar)
+    state0 = Mocca.setup_process_state(model;
+        Pressure = 1 * bar,
+        Temperature = 298.15,
+        WallTemperature = constants.T_a,
+        y = [1e-10, 1.0 - 1e-10]
+    )
+
+    parameters = Mocca.setup_process_parameters(model)
+    sim_forces, timesteps = Mocca.setup_forces(model, [t_ads], ["adsorption"];
+        num_cycles = 1, max_dt = maxdt)
+
+    case = Mocca.MoccaCase(model, timesteps, sim_forces;
+        state0 = state0, parameters = parameters)
+
+    states, timesteps_out = Mocca.simulate_process(case;
+        timestep_selector_cfg = (y = 0.01, Temperature = 10.0, Pressure = 10.0),
+        output_substates = true,
+        info_level = -1
+    )
+
+    return states, model, timesteps_out
+end
+
+function dcb_metrics(states, model, timesteps_out)
+    final = states[end]
+    nc = number_of_cells(model.domain)
+
+    breakthrough_time = nothing
+    t = 0.0
+    for (i, state) in enumerate(states)
+        t += timesteps_out[i]
+        if state[:y][1, nc] > 0.01
+            breakthrough_time = t
+            break
+        end
+    end
+
+    return (
+        y_outlet_CO2 = final[:y][1, nc],
+        y_outlet_N2 = final[:y][2, nc],
+        T_outlet = final[:Temperature][nc],
+        T_inlet = final[:Temperature][1],
+        T_max = maximum(final[:Temperature]),
+        P_outlet = final[:Pressure][nc],
+        P_inlet = final[:Pressure][1],
+        total_CO2_adsorbed = sum(final[:AdsorbedConcentration][1, :]),
+        breakthrough_time = breakthrough_time,
+    )
+end
+
+@testset "DCB Regression" begin
+    states, model, timesteps_out = run_dcb_simulation()
+    m = dcb_metrics(states, model, timesteps_out)
+
+    ref = (
+        y_outlet_CO2 = 0.1486736097381964,
+        y_outlet_N2 = 0.851326390261803,
+        T_outlet = 301.6723625454626,
+        T_inlet = 298.1499948649097,
+        T_max = 301.6723625454626,
+        P_outlet = 100004.81147504161,
+        P_inlet = 101884.19452854067,
+        total_CO2_adsorbed = 774447.598037862,
+        breakthrough_time = 484.5518204728325,
+    )
+
+    @test isapprox(m.y_outlet_CO2, ref.y_outlet_CO2, rtol=1e-3)
+    @test isapprox(m.y_outlet_N2, ref.y_outlet_N2, rtol=1e-3)
+    @test isapprox(m.T_outlet, ref.T_outlet, atol=0.5)
+    @test isapprox(m.T_inlet, ref.T_inlet, atol=0.5)
+    @test isapprox(m.T_max, ref.T_max, atol=0.5)
+    @test isapprox(m.P_outlet, ref.P_outlet, rtol=1e-3)
+    @test isapprox(m.P_inlet, ref.P_inlet, rtol=1e-3)
+    @test isapprox(m.total_CO2_adsorbed, ref.total_CO2_adsorbed, rtol=1e-3)
+    @test isapprox(m.breakthrough_time, ref.breakthrough_time, rtol=1e-3)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 668e2c5..7fda29d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -29,4 +29,9 @@ using StaticArrays
         include("jutul_integration.jl")
     end
 
+    @testset "Regression Tests" begin
+        include("regression/dcb_regression.jl")
+        include("regression/cyclic_regression.jl")
+    end
+
 end