diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000..9d75071
Binary files /dev/null and b/.DS_Store differ
diff --git a/src/.DS_Store b/src/.DS_Store
new file mode 100644
index 0000000..7b69a27
Binary files /dev/null and b/src/.DS_Store differ
diff --git a/src/pathsim_chem/.DS_Store b/src/pathsim_chem/.DS_Store
new file mode 100644
index 0000000..5b72ad7
Binary files /dev/null and b/src/pathsim_chem/.DS_Store differ
diff --git a/src/pathsim_chem/__init__.py b/src/pathsim_chem/__init__.py
index 8a72535..1f4f7c6 100644
--- a/src/pathsim_chem/__init__.py
+++ b/src/pathsim_chem/__init__.py
@@ -11,3 +11,6 @@
     __version__ = "unknown"
 
 __all__ = ["__version__"]
+
+#for direct block import from main package
+from .tritium import *
diff --git a/src/pathsim_chem/tritium/.DS_Store b/src/pathsim_chem/tritium/.DS_Store
new file mode 100644
index 0000000..5008ddf
Binary files /dev/null and b/src/pathsim_chem/tritium/.DS_Store differ
diff --git a/src/pathsim_chem/tritium/README.md b/src/pathsim_chem/tritium/README.md
new file mode 100644
index 0000000..4835a1f
--- /dev/null
+++ b/src/pathsim_chem/tritium/README.md
@@ -0,0 +1,3 @@
+# Tritium Toolbox
+
+This module contains blocks for nuclear fusion modeling. Especially for tritium cycle modeling.
\ No newline at end of file
diff --git a/src/pathsim_chem/tritium/__init__.py b/src/pathsim_chem/tritium/__init__.py
new file mode 100644
index 0000000..4cae89e
--- /dev/null
+++ b/src/pathsim_chem/tritium/__init__.py
@@ -0,0 +1,4 @@
+from .residencetime import *
+from .splitter import *
+from .bubbler import *
+# from .tcap import *
\ No newline at end of file
diff --git a/src/pathsim_chem/tritium/bubbler.py b/src/pathsim_chem/tritium/bubbler.py
new file mode 100644
index 0000000..f272858
--- /dev/null
+++ b/src/pathsim_chem/tritium/bubbler.py
@@ -0,0 +1,218 @@
+#########################################################################################
+##
+##                                      Bubbler Block
+##                                (blocks/fusion/bubbler.py)
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+
+import numpy as np
+
+from ..dynsys import DynamicalSystem
+from ...events.schedule import ScheduleList
+
+
+# BLOCK DEFIINITIONS ====================================================================
+
+class Bubbler4(DynamicalSystem):
+    """
+    Tritium bubbling system with sequential vial collection stages.
+
+    This block models a tritium collection system used in fusion reactor blanket 
+    purge gas processing. The system bubbles tritium-containing gas through a series 
+    of liquid-filled vials to capture and concentrate tritium for measurement and 
+    inventory tracking.
+
+
+    Physical Description
+    --------------------
+    The bubbler consists of two parallel processing chains:
+
+    **Soluble Chain (Vials 1-2):** 
+    Tritium already in soluble forms (HTO, HT) flows sequentially through 
+    vials 1 and 2. Each vial has a collection efficiency :math:`\\eta_{vial}`, 
+    representing the fraction of tritium that dissolves into the liquid phase 
+    and is retained.
+
+    **Insoluble Chain (Vials 3-4):** 
+    Tritium in insoluble forms (T₂, organically bound) first undergoes catalytic 
+    conversion to soluble forms with efficiency :math:`\\alpha_{conv}`. The 
+    converted tritium, along with uncaptured soluble tritium from the first chain, 
+    then flows through vials 3 and 4 with the same collection efficiency.
+
+
+    Mathematical Formulation
+    -------------------------
+    The system is governed by the following differential equations for the 
+    vial inventories :math:`x_i`:
+
+    .. math::
+
+        \\frac{dx_1}{dt} &= \\eta_{vial} \\cdot u_{sol}
+
+        \\frac{dx_2}{dt} &= \\eta_{vial} \\cdot (1-\\eta_{vial}) \\cdot u_{sol}
+
+        \\frac{dx_3}{dt} &= \\eta_{vial} \\cdot [\\alpha_{conv} \\cdot u_{insol} + (1-\\eta_{vial})^2 \\cdot u_{sol}]
+
+        \\frac{dx_4}{dt} &= \\eta_{vial} \\cdot (1-\\eta_{vial}) \\cdot [\\alpha_{conv} \\cdot u_{insol} + (1-\\eta_{vial})^2 \\cdot u_{sol}]
+
+
+    The sample output represents uncaptured tritium exiting the system:
+
+    .. math::
+    
+        y_{sample} = (1-\\alpha_{conv}) \\cdot u_{insol} + (1-\\eta_{vial})^2 \\cdot [\\alpha_{conv} \\cdot u_{insol} + (1-\\eta_{vial})^2 \\cdot u_{sol}]
+
+
+    Where:
+        - :math:`u_{sol}` = soluble tritium input flow rate
+        - :math:`u_{insol}` = insoluble tritium input flow rate  
+        - :math:`\\eta_{vial}` = vial collection efficiency
+        - :math:`\\alpha_{conv}` = conversion efficiency from insoluble to soluble
+        - :math:`x_i` = tritium inventory in vial i
+
+    Parameters
+    ----------
+    conversion_efficiency : float 
+       Conversion efficiency from insoluble to soluble forms (:math:`\\alpha_{conv}`), 
+       between 0 and 1.
+    vial_efficiency : float 
+       Collection efficiency of each vial (:math:`\\eta_{vial}`), between 0 and 1.
+    replacement_times : float | list[float] | list[list[float]] 
+        Times at which each vial is replaced with a fresh one. If None, no 
+        replacement events are created. If a single value is provided, it is 
+        used for all vials. If a single list of floats is provided, it will be 
+        used for all vials. If a list of lists is provided, each sublist 
+        corresponds to the replacement times for each vial.
+
+    Notes
+    -----
+    Vial replacement is modeled as instantaneous reset events that set the 
+    corresponding vial inventory to zero, simulating the physical replacement 
+    of a full vial with an empty one.
+    """
+
+    _port_map_out = {
+        "vial1": 0,
+        "vial2": 1,
+        "vial3": 2,
+        "vial4": 3,
+        "sample_out": 4,
+    }
+    _port_map_in = {
+        "sample_in_soluble": 0,
+        "sample_in_insoluble": 1,
+    }
+
+    def __init__(
+        self,
+        conversion_efficiency=0.9,
+        vial_efficiency=0.9,
+        replacement_times=None,
+        ):
+
+        #bubbler parameters
+        self.replacement_times = replacement_times
+        self.vial_efficiency = vial_efficiency
+        self.conversion_efficiency = conversion_efficiency
+
+        #dynamical component, ode rhs
+        def _fn_d(x, u, t):
+
+            #short
+            ve = self.vial_efficiency
+            ce = self.conversion_efficiency
+
+            #unpack inputs
+            sol, ins = u
+
+            #compute vial content change rates
+            dv1 = ve * sol
+            dv2 = dv1 * (1 - ve)
+            dv3 = ve * (ce * ins + (1 - ve)**2 * sol)
+            dv4 = dv3 * (1 - ve)
+
+            return np.array([dv1, dv2, dv3, dv4])
+
+        #algebraic output component
+        def _fn_a(x, u, t):
+
+            #short
+            ve = self.vial_efficiency
+            ce = self.conversion_efficiency
+
+            #unpack inputs
+            sol, ins = u
+
+            sample_out = (1 - ce) * ins + (1 - ve)**2 * (ce * ins + (1 - ve)**2 * sol)
+
+            return np.hstack([x, sample_out])
+
+        #initialization just like `DynamicalSystem` block
+        super().__init__(func_dyn=_fn_d, func_alg=_fn_a, initial_value=np.zeros(4))
+
+        #create internal vial reset events
+        self._create_reset_events()
+
+
+    def _create_reset_event_vial(self, i, reset_times):
+        """Define event action function and return a `ScheduleList` event 
+        per vial `i` that triggers at predefined `reset_times`. 
+        """
+
+        def reset_vial_i(_):
+            #get the full engine state
+            x = self.engine.get()
+            #set index 'i' to zero
+            x[i] = 0.0
+            #set the full engine state 
+            self.engine.set(x)
+
+        return ScheduleList(
+            times_evt=reset_times, 
+            func_act=reset_vial_i
+            )
+
+
+    def _create_reset_events(self):
+        """Create reset events for all vials based on the replacement times.
+
+        Raises
+        ------
+        ValueError : If reset_times is not valid.
+
+        Returns
+        -------
+        events : list[ScheduleList]
+            list of reset events for vials
+        """
+
+        replacement_times = self.replacement_times
+        self.events = []
+
+        # if reset_times is a single list use it for all vials
+        if replacement_times is None:
+            return 
+
+        if isinstance(replacement_times, (int, float)):
+            replacement_times = [replacement_times]
+
+        # if it's a flat list use it for all vials
+        elif isinstance(replacement_times, list) and all(
+            isinstance(t, (int, float)) for t in replacement_times
+            ):
+            replacement_times = [replacement_times] * 4
+
+        elif isinstance(replacement_times, np.ndarray) and replacement_times.ndim == 1:
+            replacement_times = [replacement_times.tolist()] * 4
+
+        elif isinstance(replacement_times, list) and len(replacement_times) != 4:
+            raise ValueError(
+                "replacement_times must be a single value or a list with the same length as the number of vials"
+            )
+
+        #create the internal events
+        self.events = [
+            self._create_reset_event_vial(i, ts) for i, ts in enumerate(replacement_times)
+            ]
diff --git a/src/pathsim_chem/tritium/residencetime.py b/src/pathsim_chem/tritium/residencetime.py
new file mode 100644
index 0000000..814c654
--- /dev/null
+++ b/src/pathsim_chem/tritium/residencetime.py
@@ -0,0 +1,123 @@
+#########################################################################################
+##
+##                           Blocks for residence time modeling
+##                            (blocks/fusion/residencetime.py)
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+
+import numpy as np
+
+from ..dynsys import DynamicalSystem
+
+
+# BLOCKS ================================================================================
+
+class ResidenceTime(DynamicalSystem):
+    """Chemical process block with residence time model.
+
+    This block implements an internal 1st order linear ode with 
+    multiple inputs, outputs, an internal constant source term 
+    and no direct passthrough.
+
+    The internal ODE with inputs :math:`u_i` :
+
+    .. math::
+        
+        \\dot{x} = - x / \\tau + \\mathrm{src} + \\sum_i \\beta_i u_i
+
+
+    And the output equation for every output `i` :
+
+    .. math::
+    
+        y_i = \\gamma_i x
+
+
+    Parameters
+    ----------
+    tau : float
+        residence time, inverse natural frequency (eigenvalue)
+    betas: None | list[float] | np.ndarray[float]
+        weights of inputs that are accumulated in state, optional
+    gammas : None | list[float] | np.ndarray[float]
+        weights of states (fractions) for output, optional
+    initial_value : float
+        initial value of state / initial quantity of process
+    source_term : float
+        constant source term / generation term of the process
+    """
+
+    def __init__(self, tau=1, betas=None, gammas=None, initial_value=0, source_term=0):
+
+        #input validation
+        if np.isclose(tau, 0):
+            raise ValueError(f"'tau' must be nonzero but is {tau}")
+    
+        #time constant and input/output weights
+        self.tau = tau
+        self.betas = 1 if betas is None else np.array(betas)
+        self.gammas = 1 if gammas is None else np.array(gammas)
+        self.source_term = source_term
+
+        #rhs of residence time ode
+        def _fn_d(x, u, t):
+            return -x/self.tau + self.source_term + sum(self.betas*u)
+
+        #jacobian of rhs wrt x
+        def _jc_d(x, u, t):
+            return -1/self.tau
+
+        #output function of residence time ode
+        def _fn_a(x, u, t):
+            return self.gammas * x
+
+        #initialization just like `DynamicalSystem` block
+        super().__init__(func_dyn=_fn_d, jac_dyn=_jc_d, func_alg=_fn_a, initial_value=initial_value)
+
+
+class Process(ResidenceTime):
+    """Simplified version of the `ResidenceTime` model block
+    with all inputs being summed equally and only the state 
+    and the flux being returned to the output
+
+    This block implements an internal 1st order linear ode with 
+    multiple inputs, outputs and no direct passthrough.
+
+    The internal ODE with inputs :math:`u_i` :
+
+    .. math::
+        
+        \\dot{x} = - x / \\tau + \\mathrm{src} + \\sum_i u_i
+
+
+    And the output equations for output `i=0` and `i=1`:
+
+    .. math::
+    
+        y_0 = x
+
+    .. math::
+
+        y_1 = x / \\tau
+
+
+    Parameters
+    ----------
+    tau : float
+        residence time, inverse natural frequency (eigenvalue)
+    initial_value : float
+        initial value of state / initial quantity of process
+    source_term : float
+        constant source term / generation term of the process
+    """
+
+    #max number of ports
+    _n_out_max = 2
+
+    #maps for input and output port labels
+    _port_map_out = {"x": 0, "x/tau": 1}
+
+    def __init__(self, tau=1, initial_value=0, source_term=0):
+        super().__init__(tau, 1, [1, 1/tau], initial_value, source_term)
diff --git a/src/pathsim_chem/tritium/splitter.py b/src/pathsim_chem/tritium/splitter.py
new file mode 100644
index 0000000..e034244
--- /dev/null
+++ b/src/pathsim_chem/tritium/splitter.py
@@ -0,0 +1,50 @@
+#########################################################################################
+##
+##                                 Flow Splitter Block 
+##                             (blocks/fusion/splitter.py)
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+
+import numpy as np
+
+from ..function import Function
+
+
+# BLOCKS ================================================================================
+
+class Splitter(Function):
+    """Splitter block that splits the input signal into multiple 
+    outputs weighted with the specified fractions.
+
+    Note
+    ----
+    The output fractions must sum to one.
+    
+    Parameters
+    ----------
+    fractions : np.ndarray | list
+        fractions to split the input signal into, 
+        must sum up to one
+    """
+
+    #max number of ports
+    _n_in_max = 1
+
+    #maps for input and output port labels
+    _port_map_in = {"in": 0}
+    
+    def __init__(self, fractions=None):
+
+        self.fractions = np.ones(1) if fractions is None else np.array(fractions)
+
+        #input validation
+        if not np.isclose(sum(self.fractions), 1):
+            raise ValueError(f"'fractions' must sum to one and not {sum(self.fractions)}")
+
+        #initialize like `Function` block
+        super().__init__(func=lambda u: self.fractions*u)
+
+        #dynamically define output port map based on fractions
+        self._port_map_out = {f"out {fr}": i for i, fr in enumerate(self.fractions)}
diff --git a/src/pathsim_chem/tritium/tcap.py b/src/pathsim_chem/tritium/tcap.py
new file mode 100644
index 0000000..8473127
--- /dev/null
+++ b/src/pathsim_chem/tritium/tcap.py
@@ -0,0 +1,27 @@
+#########################################################################################
+##
+##            Blocks for Thermal Cycle Absorption Process (TCAP) modelling
+##                              (blocks/fusion/tcap.py)
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+
+from ..ode import ODE
+
+
+# BLOCKS ================================================================================
+
+class TCAP1D(ODE):
+    """This block models the Thermal Cycle Absorption Process (TCAP) in 1d. 
+
+    The model uses a 1d finite difference spatial discretization to construct 
+    a nonlinear ODE internally as proposed in 
+
+        https://doi.org/10.1016/j.ijhydene.2023.03.101
+
+
+    """
+    raise NotImplementedError("TCAP1D block is currently not impolemented!")
+
+
diff --git a/tests/.DS_Store b/tests/.DS_Store
new file mode 100644
index 0000000..f53d21c
Binary files /dev/null and b/tests/.DS_Store differ
diff --git a/tests/tritium/__init__.py b/tests/tritium/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/tritium/test_bubbler.py b/tests/tritium/test_bubbler.py
new file mode 100644
index 0000000..8cb5163
--- /dev/null
+++ b/tests/tritium/test_bubbler.py
@@ -0,0 +1,174 @@
+########################################################################################
+##
+##                                  TESTS FOR 
+##                             'tritium.bubbler.py'
+##
+########################################################################################
+
+# IMPORTS ==============================================================================
+
+import unittest
+import numpy as np
+
+from pathsim_chem.tritium import Bubbler4
+
+from pathsim.solvers import EUF
+from pathsim.events.schedule import ScheduleList
+
+
+# TESTS ================================================================================
+
+class TestFusionBubbler4(unittest.TestCase):
+    """
+    Test the implementation of the 'Bubbler4' block class from the fusion toolbox.
+    The block inherits from `ODE` and models a 4-vial tritium collection system.
+    """
+    
+    def test_init(self):
+        """Test initialization with various parameter combinations"""
+        
+        # Default initialization
+        B = Bubbler4()
+        self.assertEqual(B.conversion_efficiency, 0.9)
+        self.assertEqual(B.vial_efficiency, 0.9)
+        self.assertEqual(B.replacement_times, None)
+        self.assertEqual(len(B.events), 0)  # No events when replacement_times is None
+        
+        # Specific initialization
+        B = Bubbler4(conversion_efficiency=0.8, vial_efficiency=0.7, replacement_times=[100, 200])
+        self.assertEqual(B.conversion_efficiency, 0.8)
+        self.assertEqual(B.vial_efficiency, 0.7)
+        self.assertEqual(B.replacement_times, [100, 200])
+        
+        # Set solver to check internal solver instance
+        B.set_solver(EUF, parent=None)
+        self.assertTrue(B.engine)
+        np.testing.assert_array_equal(B.engine.initial_value, np.zeros(4))
+    
+
+    def test_init_replacement_times(self):
+        """Test different replacement_times configurations"""
+        
+        # Single list - should replicate for all vials
+        B = Bubbler4(replacement_times=[100, 200, 300])
+        B.set_solver(EUF, parent=None)
+        self.assertEqual(len(B.events), 4)
+        
+        # List of lists - one per vial
+        times = [[100, 200], [150, 250], [120, 220], [180, 280]]
+        B = Bubbler4(replacement_times=times)
+        B.set_solver(EUF, parent=None)
+        self.assertEqual(len(B.events), 4)
+        for event in B.events:
+            self.assertIsInstance(event, ScheduleList)
+        
+        # NumPy array
+        B = Bubbler4(replacement_times=np.array([100, 200, 300]))
+        B.set_solver(EUF, parent=None)
+        self.assertEqual(len(B.events), 4)
+        
+        # Invalid case - wrong length
+        with self.assertRaises(ValueError):
+            B = Bubbler4(replacement_times=[[100], [200], [300]])  # Only 3 vials
+    
+
+    def test_update_outputs(self):
+        """Test the update method and output calculations"""
+        
+        # Test with default parameters
+        B = Bubbler4()
+        B.set_solver(EUF, parent=None)
+        
+        # Set some test inputs
+        B.inputs[0] = 10.0  # soluble input
+        B.inputs[1] = 5.0   # insoluble input
+        
+        B.update(0.0)
+        
+        # Check that outputs are set (initial state is zero)
+        self.assertEqual(B.outputs[0], 0.0)  # vial 1
+        self.assertEqual(B.outputs[1], 0.0)  # vial 2
+        self.assertEqual(B.outputs[2], 0.0)  # vial 3
+        self.assertEqual(B.outputs[3], 0.0)  # vial 4
+        
+        # Calculate expected sample_out with default parameters (ve=0.9, ce=0.9)
+        ve, ce = 0.9, 0.9
+        sol, ins = 10.0, 5.0
+        expected_sample_out = (1-ce)*ins + (1-ve)**2 * (ce*ins + (1-ve)**2 * sol)
+        self.assertAlmostEqual(B.outputs[4], expected_sample_out, places=10)
+
+
+    def test_mass_conservation(self):
+        """Test that mass is conserved in the system"""
+        
+        B = Bubbler4(conversion_efficiency=0.8, vial_efficiency=0.7)
+        B.set_solver(EUF, parent=None)
+        
+        # Set test inputs
+        sol_in, ins_in = 10.0, 5.0
+        B.inputs[0] = sol_in
+        B.inputs[1] = ins_in
+        
+        # Get rates and sample_out
+        u = np.array([sol_in, ins_in])
+        x = np.zeros(4)
+        rates = B.func_dyn(x, u, 0.0)
+        
+        B.update(0.0)
+        sample_out = B.outputs[4]
+        
+        # Total input rate
+        total_in = sol_in + ins_in
+        
+        # Total accumulation rate in vials
+        total_vial_rates = np.sum(rates)
+        
+        # Mass conservation: input = vial accumulation + sample output
+        total_out = total_vial_rates + sample_out
+        
+        self.assertAlmostEqual(total_in, total_out, places=10)
+    
+
+    def test_vial_reset_functionality(self):
+        """Test that vial reset events work correctly"""
+        
+        B = Bubbler4(replacement_times=[10.0])
+        B.set_solver(EUF, parent=None)
+        
+        # Manually set some vial inventories
+        x = np.array([1.0, 2.0, 3.0, 4.0])
+        B.engine.set(x)
+        
+        # Test reset function for vial 0
+        reset_event = B.events[0]
+        reset_func = reset_event.func_act
+        
+        # Execute reset
+        reset_func(None)
+        
+        # Check that vial 0 was reset but others remain
+        x_after = B.engine.get()
+        self.assertEqual(x_after[0], 0.0)
+        self.assertEqual(x_after[1], 2.0)
+        self.assertEqual(x_after[2], 3.0)
+        self.assertEqual(x_after[3], 4.0)
+
+        # Test reset function for vial 3
+        reset_event = B.events[3]
+        reset_func = reset_event.func_act
+        
+        # Execute reset
+        reset_func(None)
+        
+        # Check that vial 0 was reset but others remain
+        x_after = B.engine.get()
+        self.assertEqual(x_after[0], 0.0)
+        self.assertEqual(x_after[1], 2.0)
+        self.assertEqual(x_after[2], 3.0)
+        self.assertEqual(x_after[3], 0.0)
+
+
+# RUN TESTS LOCALLY ====================================================================
+
+if __name__ == '__main__':
+    unittest.main(verbosity=2)
\ No newline at end of file
diff --git a/tests/tritium/test_residencetime.py b/tests/tritium/test_residencetime.py
new file mode 100644
index 0000000..447f4ae
--- /dev/null
+++ b/tests/tritium/test_residencetime.py
@@ -0,0 +1,115 @@
+########################################################################################
+##
+##                                  TESTS FOR 
+##                          'tritium.residencetime.py'
+##
+########################################################################################
+
+# IMPORTS ==============================================================================
+
+import unittest
+import numpy as np
+
+from pathsim_chem.tritium import ResidenceTime, Process
+
+from pathsim.solvers import EUF
+
+
+# TESTS ================================================================================
+
+class TestTritiumResidenceTime(unittest.TestCase):
+    """
+    Test the implementation of the 'ResidenceTime' block class from the fusion toolbox
+
+    The block inherits from `ODE`
+    """
+    
+    def test_init(self):
+
+        #default initialization
+        R = ResidenceTime()
+
+        self.assertEqual(R.tau, 1)
+        self.assertEqual(R.betas, 1)
+        self.assertEqual(R.gammas, 1)
+        self.assertEqual(R.source_term, 0)
+
+        #input validation
+        with self.assertRaises(ValueError):
+            R = ResidenceTime(tau=0)
+
+        #specific initialization
+        R = ResidenceTime(tau=0.1, betas=[1, 3, 5], gammas=[2, 4], initial_value=10, source_term=11)
+
+        self.assertEqual(R.tau,0.1)
+        self.assertTrue(np.allclose(R.betas, np.array([1, 3, 5])))
+        self.assertTrue(np.allclose(R.gammas, np.array([2, 4])))
+
+        #set solver to check internal solver instance
+        R.set_solver(EUF, parent=None)
+
+        self.assertTrue(R.engine)
+        self.assertEqual(R.engine.initial_value, 10)
+
+
+    def test_update(self):
+
+        #default
+        R = ResidenceTime()
+        R.set_solver(EUF, parent=None)
+
+        R.update(None)
+
+        self.assertEqual(R.outputs[0], 0)
+
+        #specific
+        R = ResidenceTime(tau=0.1, betas=[1, 3, 5], gammas=[2, 4], initial_value=11, source_term=12)
+        R.set_solver(EUF, parent=None)
+
+        R.update(None)
+
+        self.assertEqual(R.outputs[0], 22) # initial_value * gammas[0]
+        self.assertEqual(R.outputs[1], 44) # initial_value * gammas[1]
+
+
+
+class TestTritiumProcess(unittest.TestCase):
+    """
+    Test the implementation of the 'Process' block class from the fusion toolbox.
+
+    This block inherits from `ResidenceTime`, just testing the initialization.
+    """
+
+    def test_init(self):
+
+        #default initialization
+        P = Process()
+
+        self.assertEqual(P.tau, 1)
+        self.assertEqual(P.betas, 1)
+        self.assertTrue(np.allclose(P.gammas, [1, 1])) # 1, 1/tau
+        self.assertEqual(P.source_term, 0)
+
+        #specific initialization
+        P = Process(tau=0.1, source_term=33, initial_value=44)
+
+        self.assertEqual(P.tau,0.1)
+        self.assertEqual(P.betas, 1)
+        self.assertTrue(np.allclose(P.gammas, np.array([1, 10])))
+        self.assertEqual(P.source_term, 33)
+
+        #set solver to check internal solver instance
+        P.set_solver(EUF, parent=None)
+
+        self.assertTrue(P.engine)
+        self.assertEqual(P.engine.initial_value, 44)
+
+
+
+
+
+
+# RUN TESTS LOCALLY ====================================================================
+
+if __name__ == '__main__':
+    unittest.main(verbosity=2)
diff --git a/tests/tritium/test_splitter.py b/tests/tritium/test_splitter.py
new file mode 100644
index 0000000..71f3e97
--- /dev/null
+++ b/tests/tritium/test_splitter.py
@@ -0,0 +1,61 @@
+########################################################################################
+##
+##                                  TESTS FOR 
+##                            'tritium.splitter.py'
+##
+########################################################################################
+
+# IMPORTS ==============================================================================
+
+import unittest
+import numpy as np
+
+from pathsim_chem.tritium import Splitter
+
+
+# TESTS ================================================================================
+
+class TestTritiumSplitter(unittest.TestCase):
+    """
+    Test the implementation of the 'Splitter' block class from the fusion toolbox
+    """
+
+    def test_init(self):
+        
+        #default initialization
+        S = Splitter()
+        self.assertEqual(S.fractions, np.ones(1))
+
+        #input validation
+        for fracs in [[1, 3], [0.4, 0.6, 0.001], [0.33, 0.33, 0.33]]:
+            with self.assertRaises(ValueError):
+                S = Splitter(fracs)
+
+        #special initialization
+        S = Splitter([0.4, 0.5, 0.1])
+        self.assertEqual(sum(S.fractions - np.array([0.4, 0.5, 0.1])), 0)
+
+        #test the automatic port maps
+        self.assertEqual(S._port_map_out, {"out 0.4":0, "out 0.5":1, "out 0.1":2})
+
+
+    def test_update(self):
+
+        S = Splitter([0.4, 0.5, 0.1])
+
+        #set block inputs
+        S.inputs[0] = 2
+
+        #update block
+        S.update(None)
+
+        #test if update was correct
+        self.assertEqual(S.outputs[0], 0.8)
+        self.assertEqual(S.outputs[1], 1)
+        self.assertEqual(S.outputs[2], 0.2)
+
+
+# RUN TESTS LOCALLY ====================================================================
+
+if __name__ == '__main__':
+    unittest.main(verbosity=2)