diff --git a/book/_toc.yml b/book/_toc.yml
index 0ae07d8..ee441bb 100644
--- a/book/_toc.yml
+++ b/book/_toc.yml
@@ -21,6 +21,13 @@ parts:
     - file: content/material/material_htm
     - file: content/material/material_advanced
 
+  - caption: Boundary Conditions
+    chapters: 
+    - file: content/boundary_conditions/concentration
+    - file: content/boundary_conditions/fluxes
+    - file: content/boundary_conditions/surface_reactions
+    - file: content/boundary_conditions/heat_transfer
+
   - caption: Species & Reactions
     chapters:
     - file: content/species_reactions/species
diff --git a/book/content/boundary_conditions/concentration.md b/book/content/boundary_conditions/concentration.md
new file mode 100644
index 0000000..ec09ff6
--- /dev/null
+++ b/book/content/boundary_conditions/concentration.md
@@ -0,0 +1,169 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.18.1
+kernelspec:
+  display_name: festim-workshop
+  language: python
+  name: python3
+---
+
+# Concentration #
+
+Boundary conditions (BCs) are essential to FESTIM simulations, as they describe the mathematical problem at the boundaries of the simulated domain. Read more about BCs _[here](https://festim.readthedocs.io/en/fenicsx/userguide/boundary_conditions.html)_.
+
+This tutorial goes over how to define concentration boundary conditions for hydrogen transport simulations.
+
+Objectives:
+* Define a fixed concentration BC
+* Choose which solubility law (Sieverts' or Henry's)
+* Solve a hydrogen transport problem with plasma implantation
+
++++
+
+## Defining fixed concentration ##
+
+BCs need to be assigned to surfaces using FESTIM's `SurfaceSubdomain` class. To define the concentration of a specific species, we can use `FixedConcentrationBC`:
+
+```{code-cell} ipython3
+from festim import FixedConcentrationBC, Species, SurfaceSubdomain
+
+boundary = SurfaceSubdomain(id=1)
+H = Species(name="Hydrogen")
+
+my_bc = FixedConcentrationBC(subdomain=boundary, value=10, species=H)
+```
+
+The imposed concentration can be dependent on space, time and temperature, such as:
+
+$$ 
+
+BC = 10 + x^2 + T(t+2)
+
+$$
+
+```{code-cell} ipython3
+my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2)
+
+my_bc = FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H)
+```
+
+## Choosing a solubility law ##
+
+Users can define the surface concentration using either Sieverts’ law, $c = S(T)\sqrt P$, or Henry's law, $c=K_H P$, where $S(T)$ and $K_H$ denote the temperature-dependent Sieverts’ and Henry’s solubility coefficients, respectively, and $P$ is the partial pressure of the species on the surface. 
+
+For Sieverts' law of solubility, we can use `festim.SievertsBC`:
+
+```{code-cell} ipython3
+from festim import SievertsBC, SurfaceSubdomain, Species
+
+boundary = SurfaceSubdomain(id=1)
+H = Species(name="Hydrogen")
+
+custom_pressure_value = lambda t: 2 + t
+my_bc = SievertsBC(subdomain=3, S_0=2, E_S=0.1, species=H, pressure=custom_pressure_value)
+```
+
+Similarly, for Henry's law of solubility, we can use `festim.HenrysBC`:
+
+```{code-cell} ipython3
+from festim import HenrysBC
+
+pressure_value = lambda t: 5 * t
+my_bc = HenrysBC(subdomain=3, H_0=1.5, E_H=0.2, species=H, pressure=pressure_value)
+```
+
+## Plasma implantation approximation ##
+
++++
+
+We can also approximate plasma implantation using FESTIM's `ParticleSource` class, which is helpful in modeling thermo-desorption spectra (TDS) experiments or including the effect of plasma exposure on hydrogen transport. Learn more about how FESTIM approximates plasma implantation _[here](https://festim.readthedocs.io/en/fenicsx/theory.html)_.
+
+Consider the following 1D plasma implantation problem, where we represent the plasma as a hydrogen source $S_{ext}$:
+
+$$ S_{ext} = \varphi \cdot f(x) $$
+
+$$\varphi = 1\cdot 10^{13} \quad \mathrm{m}^{-2}\mathrm{s}^{-1}$$
+
+where  $\varphi$ is the implantation flux and $f(x)$ is a Gaussian spatial distribution (distribution mean value of 0.5 $\text{m}$ and width of 1 $\text{m}$).
+
+First, we setup a 1D mesh ranging from $ [0,1] $ and assign the subdomains and material:
+
+```{code-cell} ipython3
+import festim as F
+import ufl
+import numpy as np
+
+my_model = F.HydrogenTransportProblem()
+vertices = np.linspace(0,1,2000)
+my_model.mesh = F.Mesh1D(vertices)
+
+mat = F.Material(D_0=0.1, E_D=0.01)
+
+volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat)
+left_boundary = F.SurfaceSubdomain1D(id=1, x=0)
+right_boundary = F.SurfaceSubdomain1D(id=2, x=1)
+
+my_model.subdomains = [
+    volume_subdomain,
+    left_boundary,
+    right_boundary,
+]
+```
+
+Now, we define our `incident_flux` and `gaussian_distribution` function. We can use `ParticleSource` to represent the source term, and then add it to our model:
+
+```{code-cell} ipython3
+incident_flux = 1e13 
+
+def gaussian_distribution(x, center, width):
+    return (
+        1
+        / (width * (2 * ufl.pi) ** 0.5)
+        * ufl.exp(-0.5 * ((x[0] - center) / width) ** 2)
+    )
+H = F.Species("H")
+my_model.species = [H]
+
+source_term = F.ParticleSource(
+    value=lambda x: incident_flux * gaussian_distribution(x, .5, 1),
+    volume=volume_subdomain,
+    species=H,
+)
+
+my_model.sources = [source_term]
+```
+
+Finally, we assign boundary conditions (zero concentration at $x=0$ and $x=1$) and solve our problem:
+
+```{code-cell} ipython3
+my_model.boundary_conditions = [
+    F.FixedConcentrationBC(subdomain=left_boundary, value=0, species=H),
+    F.FixedConcentrationBC(subdomain=right_boundary, value=0, species=H),
+]
+
+my_model.temperature = 400
+my_model.settings = F.Settings(atol=1e10, rtol=1e-10, transient=False)
+
+profile_export = F.Profile1DExport(field=H,subdomain=volume_subdomain)
+my_model.exports = [profile_export]
+
+my_model.initialise()
+my_model.run()
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import matplotlib.pyplot as plt
+
+x = my_model.exports[0].x
+c = my_model.exports[0].data[0][0]
+
+plt.plot(x, c)
+plt.show()
+```
diff --git a/book/content/boundary_conditions/fluxes.md b/book/content/boundary_conditions/fluxes.md
new file mode 100644
index 0000000..493dedd
--- /dev/null
+++ b/book/content/boundary_conditions/fluxes.md
@@ -0,0 +1,195 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.18.1
+kernelspec:
+  display_name: festim-workshop
+  language: python
+  name: python3
+---
+
+# Particle flux #
+
+This tutorial goes over how to define particle flux boundary conditions in FESTIM.
+
+Objectives:
+* Learn how to define fluxes at boundaries
+* Impose boundary fluxes as a function of the concentration
+* Solve a problem with multi-species flux boundary conditions
+
++++
+
+## Defining flux at boundaries ##
+
+Users can impose the particle flux (Neumann BC) at boundaries using `ParticleFluxBC` class:
+
+```{code-cell} ipython3
+from festim import ParticleFluxBC, Species, SurfaceSubdomain
+
+boundary = SurfaceSubdomain(id=1)
+H = Species(name="Hydrogen")
+
+my_flux_bc = ParticleFluxBC(subdomain=boundary, value=2, species=H)
+```
+
+## Defining concentration-dependant fluxes ##
+
+Similar to the concentration, the flux can be dependent on space, time and temperature. But for particle fluxes, the values can also be dependent on a species’ concentration. 
+
+For example, let's define a hydrogen flux `J` that depends on the hydrogen concentration `c` and time `t`:
+
+$$ J(c,t) = 10t^2 + 2c $$
+
+```{code-cell} ipython3
+from festim import ParticleFluxBC, Species, SurfaceSubdomain
+
+boundary = SurfaceSubdomain(id=1)
+H = Species(name="Hydrogen")
+
+J = lambda t, c: 10*t**2 + 2*c
+
+my_flux_bc = ParticleFluxBC(
+    subdomain=boundary,
+    value=J,
+    species=H,
+    species_dependent_value={"c": H},
+)
+```
+
+## Multi-species flux boundary conditions ##
+
++++
+
+Users may also need to impose a flux boundary condition which depends on the concentrations of multiple species. 
+
+Consider the following 1D example with three species, $\mathrm{A}$, $\mathrm{B}$, and $\mathrm{C}$, with recombination fluxes $\phi_{AB}$ and $\phi_{ABC}$ that depend on the concentrations $\mathrm{c}$:
+
+$$ \phi_{AB} = -c_A c_B $$
+
+$$ \phi_{ABC}= -10 c_A c_B c_C$$
+
+We must first define each species using `Species` and then create the dictionary to be passed into `species_dependent_value`. The dictionary maps each argument in the custom flux function to the corresponding defined species. We also define our custom functions for the fluxes:
+
+```{code-cell} ipython3
+import festim as F
+
+my_model = F.HydrogenTransportProblem()
+
+A = F.Species(name="A")
+B = F.Species(name="B")
+C = F.Species(name="C")
+
+my_model.species = [A, B, C]
+species_dependent_value = {"c_A": A, "c_B": B, "c_C": C}
+
+recombination_flux_AB = lambda c_A, c_B, c_C: -c_A*c_B
+recombination_flux_ABC = lambda c_A, c_B, c_C: -10*c_A*c_B*c_C
+```
+
+Now, we create our 1D mesh and assign boundary conditions (recombination on the right, fixed concentration on the left).
+
+The boundary condition `ParticleFluxBC` must be added for each species:
+
+```{code-cell} ipython3
+import numpy as np
+
+my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100))
+
+D = 1e-2
+mat = F.Material(
+    D_0={A: 8*D, B: 7*D, C: D},
+    E_D={A: 0.01, B: 0.01, C: 0.01},
+)
+
+bulk = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat)
+left = F.SurfaceSubdomain1D(id=1, x=0)
+right = F.SurfaceSubdomain1D(id=2, x=1)
+my_model.subdomains = [bulk, left, right]
+
+my_model.boundary_conditions = [
+    F.ParticleFluxBC(
+        subdomain=right,
+        value=recombination_flux_AB,
+        species=A,
+        species_dependent_value=species_dependent_value,
+    ),
+    F.ParticleFluxBC(
+        subdomain=right,
+        value=recombination_flux_AB,
+        species=B,
+        species_dependent_value=species_dependent_value,
+    ),
+    F.ParticleFluxBC(
+        subdomain=right,
+        value=recombination_flux_ABC,
+        species=A,
+        species_dependent_value=species_dependent_value,
+    ),
+    F.ParticleFluxBC(
+        subdomain=right,
+        value=recombination_flux_ABC,
+        species=B,
+        species_dependent_value=species_dependent_value,
+    ),
+    F.ParticleFluxBC(
+        subdomain=right,
+        value=recombination_flux_ABC,
+        species=C,
+        species_dependent_value=species_dependent_value,
+    ),
+    F.FixedConcentrationBC(subdomain=left,value=1,species=A),
+    F.FixedConcentrationBC(subdomain=left,value=1,species=B),
+    F.FixedConcentrationBC(subdomain=left,value=1,species=C),     
+]
+
+right_flux_A = F.SurfaceFlux(field=A,surface=right)
+right_flux_B = F.SurfaceFlux(field=B,surface=right)
+right_flux_C = F.SurfaceFlux(field=C,surface=right)
+
+my_model.exports = [
+    right_flux_A,
+    right_flux_B,
+    right_flux_C,
+]
+```
+
+```{note}
+The diffusivity pre-factor `D_0` and activation energy `E_D` must be defined for each species in `Material`. Learn more about defining multi-species material properties __[here](https://festim-workshop.readthedocs.io/en/festim2/content/material/material_advanced.html#defining-species-dependent-material-properties)__. 
+```
+
++++
+
+Finally, let's solve and plot the profile for each species:
+
+```{code-cell} ipython3
+my_model.temperature = 300
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
+
+my_model.initialise()
+my_model.run()
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import matplotlib.pyplot as plt
+
+def plot_profile(species, **kwargs):
+    c = species.post_processing_solution.x.array[:]
+    x = species.post_processing_solution.function_space.mesh.geometry.x[:,0]
+    return plt.plot(x, c, **kwargs)
+
+for species in my_model.species:
+    plot_profile(species, label=species.name)
+
+plt.xlabel('Position')
+plt.ylabel('Concentration')
+plt.legend()
+plt.show()
+```
+
+We see the higher recombination flux for $\mathrm{ABC}$ decreases the concentration of $\mathrm{C}$ throughout the material.
diff --git a/book/content/boundary_conditions/heat_transfer.md b/book/content/boundary_conditions/heat_transfer.md
new file mode 100644
index 0000000..e78df30
--- /dev/null
+++ b/book/content/boundary_conditions/heat_transfer.md
@@ -0,0 +1,64 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.18.1
+---
+
+# Heat transfer #
+
+This tutorial goes over how to define boundary conditions for heat transfer simulations.
+
+Objectives:
+* Define homogenous temperature boundary conditions
+* Define heat flux boundary conditions
+
++++
+
+## Imposing homogenous temperature boundary conditions ##
+
+The temperature can be imposed on boundaries using `FixedTemperatureBC`:
+
+```{code-cell}
+from festim import FixedTemperatureBC, SurfaceSubdomain
+
+boundary = SurfaceSubdomain(id=1)
+my_bc = FixedTemperatureBC(subdomain=boundary, value=10)
+```
+
+To define the temperature as space or time dependent, a function can be passed to the `value` argument, such as:
+
+$$ \text{BC} = 10 + x^2 + t $$
+
+```{code-cell}
+from festim import FixedTemperatureBC, SurfaceSubdomain
+
+boundary = SurfaceSubdomain(id=1)
+BC = lambda x, t: 10 + x[0]**2 + t
+
+my_bc = FixedTemperatureBC(subdomain=boundary, value=BC)
+```
+
+## Imposing heat flux boundary conditions ##
+
++++
+
+Heat fluxes can be imposed on boundaries using `HeatFluxBC`, which can depend on space, time, and temperature, such as:
+
+$$ \text{BC} = 2x + 10t + T $$
+
+```{code-cell}
+from festim import HeatFluxBC, SurfaceSubdomain
+
+boundary = SurfaceSubdomain(id=1)
+BC = lambda x, t, T: 2 * x[0] + 10 * t + T
+
+my_flux_bc = HeatFluxBC(subdomain=boundary, value=BC)
+```
+
+```{note}
+Read more about heat transfer settings __[here](https://festim-workshop.readthedocs.io/en/festim2/content/temperatures/temperatures_advanced)__.
+```
diff --git a/book/content/boundary_conditions/surface_reactions.md b/book/content/boundary_conditions/surface_reactions.md
new file mode 100644
index 0000000..8764565
--- /dev/null
+++ b/book/content/boundary_conditions/surface_reactions.md
@@ -0,0 +1,337 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.18.1
+kernelspec:
+  display_name: festim-workshop
+  language: python
+  name: python3
+---
+
+# Surface reactions #
+
++++
+
+FESTIM V2.0 allows users to impose surface reactions on boundaries. Learn more about surface reactions __[here](https://festim-workshop.readthedocs.io/en/festim2/content/species_reactions/surface_reactions.html)__.
+
+Objectives:
+* Defining a simple surface reaction boundary condition 
+* Recombination and dissociation
+* Isotopic exhange
+* Complex isotopic exchange with multple hydrogenic species
+
++++
+
+## Defining a simple surface reaction boundary condition ##
+
+Surface reactions between adsorbed species and gas-phase products can be represented using the `SurfaceReactionBC` class. This boundary condition defines a reversible reaction of the form:
+
+$$
+A + B \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad C
+$$
+
+where $\mathrm{A}$ and $\mathrm{B}$ are surface reactants and $\mathrm{C}$ is the product species in the gas phase.
+
+### Reaction Rate Formulation ###
+
+The forward and backward rate constants follow Arrhenius laws:
+
+$$
+K_r = k_{r0} e^{-E_{kr} / (k_B T)}, \qquad
+K_d = k_{d0} e^{-E_{kd} / (k_B T)}
+$$
+
+The net surface reaction rate is given by:
+
+$$
+K = K_r c_A c_B - K_d P_C
+$$
+
+where: 
+- $\mathrm{k_{r0}}, \mathrm{k_{d0}}$ are rate pre-exponential constants  
+- $\mathrm{c_A}, \mathrm{c_B}$ are concentrations of reactant species $\mathrm{A}$ and $\mathrm{B}$ at the surface  
+- $\mathrm{P_C}$ is the partial pressure of the product species $\mathrm{C}$  
+- $\mathrm{k_B}$ is the Boltzmann constant  
+- $\mathrm{T}$ is the surface temperature  
+
+The flux of species $\mathrm{A}$ entering the surface is equal to $\mathrm{K}$ (if $\mathrm{A=B}$, the total particle flux entering the surface becomes $\mathrm{2K}$).
+
+We can use the `SurfaceReactionBC` class to impose the surface reaction above by specifying the reactants (`reactant`), gas pressure (`gas_pressure`), forward/backward rate constants (`k_r0` and `k_d0`), and rate activation energies (`E_kr` and `E_kd`):
+
+```{code-cell} ipython3
+from festim import Species, SurfaceReactionBC, SurfaceSubdomain
+
+boundary = SurfaceSubdomain(id=1)
+A = Species("A")
+B = Species("B")
+C = Species("C")
+
+my_bc = SurfaceReactionBC(
+    reactant=[A, B],
+    gas_pressure=1e5,
+    k_r0=1,
+    E_kr=0.1,
+    k_d0=0,
+    E_kd=0,
+    subdomain=boundary,
+)
+```
+
+## Recombination and dissociation ##
+
+Hydrogen recombination/dissociation can also be modelled using `SurfaceReactionBC`, such as this reaction:
+
+$$ \mathrm{H + H} \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad \mathrm{H_2} $$
+
+```{code-cell} ipython3
+from festim import Species, SurfaceReactionBC, SurfaceSubdomain
+
+boundary = SurfaceSubdomain(id=1)
+H = Species("H")
+
+my_bc = SurfaceReactionBC(
+    reactant=[H, H],
+    gas_pressure=1e5,
+    k_r0=1,
+    E_kr=0.1,
+    k_d0=1e-5,
+    E_kd=0.1,
+    subdomain=boundary,
+)
+```
+
+## Isotopic exchange ##
+
+Isotopic exchange can occur where isotopes can swap positions with other isotopes, such as:
+
+$$\mathrm{T + H_2} \rightleftharpoons \mathrm{H} + \mathrm{HT}$$
+
+These exchange processes occur through complex mechanisms and are important in understanding the behavior of hydrogenic species:
+
+If the $\mathrm{H_2}$ concentration is assumed much larger than $\mathrm{HT}$, the flux $\phi_T$ reduces to a first-order process in $\mathrm{T}$:
+
+$$
+\phi_T = 
+-\mathrm{K_{r0}} \exp\!\left(-\frac{\mathrm{E_{Kr}}}{\mathrm{k_B T}}\right) c_T^2
+-\mathrm{K_{r0}^\ast} \exp\!\left(-\frac{\mathrm{E_{Kr}^\ast}}{\mathrm{k_B T}}\right) c_{H_2} c_T
+$$
+
+Such fluxes can be implemented using `festim.ParticleFluxBC` with user-defined expressions.
+
+```{code-cell} ipython3
+import ufl
+import festim as F
+
+Kr_0 = 1.0
+E_Kr = 0.1
+
+def my_custom_recombination_flux(c, T):
+    Kr_0_custom = 1.0
+    E_Kr_custom = 0.5  # eV
+    h2_conc = 1e25  # assumed constant H2 concentration in
+
+    recombination_flux = (
+        -(Kr_0 * ufl.exp(-E_Kr / (F.k_B * T))) * c**2
+        - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * h2_conc * c
+    )
+    return recombination_flux
+```
+
+```{note}
+We can also use `F.SurfaceReactionBC` to define recombination fluxes, which we do in the next section for more complex isotopic exchanges.
+```
+
++++
+
+Let's work through a 1D example to illustrate the effect of isotopic exchange. We also define a fixed concentration on the left surface:
+
+```{code-cell} ipython3
+import numpy as np
+
+my_model = F.HydrogenTransportProblem()
+my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 100))
+
+left_surf = F.SurfaceSubdomain1D(id=1, x=0)
+right_surf = F.SurfaceSubdomain1D(id=2, x=1)
+
+material = F.Material(D_0=1e-2, E_D=0)
+
+vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material)
+
+my_model.subdomains = [vol, left_surf, right_surf]
+tritium = F.Species("T")
+my_model.species = [tritium]
+
+my_custom_flux = F.ParticleFluxBC(
+    value=my_custom_recombination_flux,
+    subdomain=right_surf,
+    species_dependent_value={"c": tritium},
+    species=tritium,
+)
+
+my_model.boundary_conditions = [
+    my_custom_flux,
+    F.FixedConcentrationBC(subdomain=left_surf, value=1, species=tritium),
+]
+
+my_model.temperature = 300
+right_flux = F.SurfaceFlux(field=tritium, surface=right_surf)
+
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
+
+my_model.initialise()
+my_model.run()
+
+data_with_H2 = (tritium.solution.x.array)
+```
+
+We can compare the surface flux to the case where we have no isotopic exchange by removing the custom boundary condition:
+
+```{code-cell} ipython3
+my_model.boundary_conditions = [
+    F.FixedConcentrationBC(subdomain=left_surf, value=1, species=tritium),
+]
+
+my_model.initialise()
+my_model.run()
+data_no_H2 = (tritium.solution.x.array)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import matplotlib.pyplot as plt
+
+x = my_model.mesh.mesh.geometry.x[:, 0]
+plt.plot(x, data_no_H2, label="No H_2")
+plt.plot(x, data_with_H2, label="With H_2")
+
+plt.xlabel("x (m)")
+plt.ylabel("Surface Flux (right)")
+plt.legend()
+plt.show()
+```
+
+We see that diffusion is present with isotopic exchange!
+
++++
+
+## Complex isotopic exchange with multple hydrogenic species ##
+
+Surface reactions can involve multiple hydrogen isotopes, allowing for the modeling of complex isotope-exchange mechanisms between species. For example, in a system with both mobile hydrogen and deuteriun, various molecular recombination pathways may occur at the surface, resulting in the formation of $H_2$, $D_2$, and $HD$:
+
+$$ \text{Reaction 1}: \mathrm{H+D} \rightarrow \mathrm{HD} \longrightarrow \phi_1 = K_{r1} c_H c_D - K_{d1}P_{HD} $$
+$$ \text{Reaction 2}: \mathrm{D+D} \rightarrow \mathrm{D_2} \longrightarrow \phi_2 = 2K_{r2} c_D^2 - K_{d2}P_{D2} $$
+$$ \text{Reaction 3}: \mathrm{D+H_2} \rightarrow \mathrm{HD + H} \longrightarrow \phi_3 = K_{r3} c_H c_D - K_{d3}P_{HD} $$
+$$ \text{Reaction 4}: \mathrm{H+H} \rightarrow \mathrm{H_2} \longrightarrow \phi_4 = 2K_{r4} c_H^2 - K_{d4}P_{H2} $$
+
+Now consider the case where deuterium diffuses from left to right and reacts with background 
+$\mathrm{H_2}$, while $\mathrm{P_{HD}}$ and $\mathrm{P_{D_2}}$ are negligible at the surface. 
+Formation of $\mathrm{H}$ at the right boundary induces back-diffusion toward the left, 
+even though none existed initially. 
+
+The boundary conditions for this scenario are:
+
+$$
+c_D(x=0) = 1, \qquad c_H(x=0) = 0, \qquad P_{H2}(x=1) = \text{1000 Pa}
+$$
+
+First, let's define a 1D mesh ranging from $\mathrm{x=[0,1]}$:
+
+```{code-cell} ipython3
+import numpy as np
+import festim as F
+
+my_model = F.HydrogenTransportProblem()
+my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 100))
+
+left_surf = F.SurfaceSubdomain1D(id=1, x=0)
+right_surf = F.SurfaceSubdomain1D(id=2, x=1)
+
+material = F.Material(D_0=1e-2, E_D=0)
+vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material)
+
+my_model.subdomains = [vol, left_surf, right_surf]
+```
+
+Now, we define our species at recombination reactions using `SurfaceReactionBC`:
+
+```{code-cell} ipython3
+H = F.Species("H")
+D = F.Species("D")
+my_model.species = [H, D]
+
+H2 = F.SurfaceReactionBC(
+    reactant=[H, H],
+    gas_pressure=100000,
+    k_r0=1,
+    E_kr=0.1,
+    k_d0=1e-5,
+    E_kd=0.1,
+    subdomain=right_surf,
+)
+
+HD = F.SurfaceReactionBC(
+    reactant=[H, D],
+    gas_pressure=0,
+    k_r0=1,
+    E_kr=0.1,
+    k_d0=1e-5,
+    E_kd=0.1,
+    subdomain=right_surf,
+)
+
+D2 = F.SurfaceReactionBC(
+    reactant=[D, D],
+    gas_pressure=0,
+    k_r0=1,
+    E_kr=0.1,
+    k_d0=1e-5,
+    E_kd=0.1,
+    subdomain=right_surf,
+)
+```
+
+Finally, we add our boundary conditions and solve the steady-state problem:
+
+```{code-cell} ipython3
+my_model.boundary_conditions = [
+    H2,
+    D2,
+    HD,
+    F.FixedConcentrationBC(subdomain=left_surf, value=1, species=D),
+    F.FixedConcentrationBC(subdomain=left_surf, value=0, species=H),
+]
+
+my_model.temperature = 300
+
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
+
+my_model.initialise()
+my_model.run()
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import matplotlib.pyplot as plt
+
+def plot_profile(species, **kwargs):
+    c = species.post_processing_solution.x.array[:]
+    x = species.post_processing_solution.function_space.mesh.geometry.x[:,0]
+    return plt.plot(x, c, **kwargs)
+
+for species in my_model.species:
+    plot_profile(species, label=species.name)
+
+plt.xlabel('Position')
+plt.ylabel('Concentration')
+plt.legend()
+plt.show()
+```
+
+We see that the background $\mathrm{H_2}$ reacts with the $\mathrm{D}$, removing the total amount of $\mathrm{D}$ from the surface. Conversely, the $\mathrm{H}$ diffuses from the surface towards the left.