diff --git a/book/_toc.yml b/book/_toc.yml
index ef38ce6..29efa33 100644
--- a/book/_toc.yml
+++ b/book/_toc.yml
@@ -21,6 +21,14 @@ parts:
     - file: content/material/material_htm
     - file: content/material/material_advanced
 
+  - caption: Boundary Conditions
+    chapters: 
+    - file: content/boundary_conditions/bc_intro
+    - file: content/boundary_conditions/h_transport_basic
+    - file: content/boundary_conditions/h_transport_advanced
+    - file: content/boundary_conditions/examples
+    - file: content/boundary_conditions/heat_transfer
+
   - caption: Species & Reactions
     chapters:
     - file: content/species_reactions/species
diff --git a/book/content/boundary_conditions/bc_intro.md b/book/content/boundary_conditions/bc_intro.md
new file mode 100644
index 0000000..b21ac27
--- /dev/null
+++ b/book/content/boundary_conditions/bc_intro.md
@@ -0,0 +1,79 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.19.1
+---
+
+# Fixed value and flux boundary condition mathematics
+
+This section discusses the math behind fixed temperature/concentration and flux boundary conditions (BCs).
+
++++
+
+## Fixed concentration/temperature boundary conditions
+
+A fixed concentration (Dirichlet) boundary condition prescribes the value of the mobile hydrogen isotope concentration at a boundary. This enforces the concentration to remain constant in time and space on the specified boundary, independent of the solution in the bulk. This also applies to a fixed temperature boundary condition.
+
+For hydrogen transport, this boundary condition is typically used to represent surfaces in equilibrium with a gas, imposed implantation conditions, or experimentally controlled concentrations. Fixed temperatures are commonly used to model infinite heat reservoirs or sinks.
+
+### Mathematical formulation
+
+On a boundary $\Gamma_D$, the mobile concentration satisfies
+
+$$
+c(\mathbf{x}, t) = c_0 \quad \text{for } \mathbf{x} \in \Gamma_D,
+$$
+
+where $c_0$ is the prescribed concentration value.
+
+This condition is enforced by directly constraining the degrees of freedom associated with the boundary.
+
+## Flux boundary conditions
+
+A particle flux (Neumann) boundary condition prescribes the normal flux of mobile hydrogen isotopes across a boundary. Unlike fixed concentration conditions, flux boundary conditions do not directly constrain the concentration value at the surface. Instead, they control the rate at which particles enter or leave the domain.
+
+Flux boundary conditions are commonly used to represent implantation from a plasma, outgassing to vacuum, permeation through a surface, or symmetry boundaries where no net transport occurs.
+
+This also applies for heat flux boundary conditions. 
+
+### Mathematical formulation
+
+The particle flux $\mathbf{J}$ is typically given by Fick’s law,
+
+$$
+\mathbf{J} = -D \nabla c,
+$$
+
+where $D$ is the diffusion coefficient and $c$ is the mobile concentration.
+
+On a boundary $\Gamma_N$, a prescribed normal flux is imposed as
+
+$$
+\mathbf{J} \cdot \mathbf{n} = g(\mathbf{x}, t) \quad \text{for } \mathbf{x} \in \Gamma_N,
+$$
+
+where:
+- $\mathbf{n}$ is the outward unit normal vector,
+- $g(\mathbf{x}, t)$ is the imposed particle flux (positive for outward flux, negative for inward flux).
+
+A special case is the **zero-flux (symmetry) boundary condition**, given by
+
+$$
+\mathbf{J} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma,
+$$
+
+which implies no net particle transport across the boundary.
+
+#### Weak formulation
+
+In the weak form, flux boundary conditions appear naturally as surface integrals after integration by parts. For a test function $v$, the boundary contribution is
+
+$$
+\int_{\Gamma_N} g(\mathbf{x}, t)\, v \, \mathrm{d}\Gamma.
+$$
+
+Because of this, flux boundary conditions are often referred to as *natural boundary conditions* and do not require explicit modification of the solution space, unlike Dirichlet conditions.
diff --git a/book/content/boundary_conditions/examples.md b/book/content/boundary_conditions/examples.md
new file mode 100644
index 0000000..979631f
--- /dev/null
+++ b/book/content/boundary_conditions/examples.md
@@ -0,0 +1,361 @@
+---
+jupytext:
+  formats: md:myst,ipynb
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.20.0
+kernelspec:
+  display_name: festim-workshop
+  language: python
+  name: python3
+---
+
+# Hydrogen transport: examples
+
+This section discusses a few advanced examples using various boundary conditions to help users understand how BCs can help model their systems.
+
+Objective:
+* Model plasma implantation using volumetric sources and approximations
+* Model complex multi-species isotopic exchange surface reactions
+
++++
+
+## Modeling plasma implantation
+
+We can model plasma implantation using FESTIM's `ParticleSource` class, which is a class used to define volumetric source terms. This is helpful in modeling thermo-desorption spectra (TDS) experiments or including the effect of plasma exposure on hydrogen transport. 
+
+Consider the following 1D plasma implantation problem, where we represent the plasma as a hydrogen source $S_{ext}$ implanted on a tungsten material that is 20 microns thick:
+
+$$ S_{ext} = \varphi_{imp} \cdot f(x) $$
+$$ L = 20 \mu \mathrm{m}$$
+$$\varphi_{imp} = 10^{13} \quad \mathrm{m}^{-2}\mathrm{s}^{-1}$$
+
+where  $\varphi_{imp}$ is the implantation flux and $f(x)$ is a Gaussian spatial distribution (distribution mean value represents implantation depth).
+
+First, we setup a 1D mesh ranging from $ [0,L] $ and assign the subdomains and material properties for tungsten:
+
+```{code-cell} ipython3
+import festim as F
+import ufl
+import numpy as np
+
+L = 20e-6
+my_model = F.HydrogenTransportProblem()
+
+vertices = np.concatenate(
+    [
+        np.linspace(0, 30e-9, num=200),
+        np.linspace(3e-6, L, num=500),
+    ]
+)
+
+my_model.mesh = F.Mesh1D(vertices)
+```
+
+```{note}
+We break up the mesh region into two regions so we can refine the region closer to the implantation depth (defined below)
+```
+
+```{code-cell} ipython3
+tungsten = F.Material(
+    D_0=4.1e-07,  # m2/s
+    E_D=0.39,  # eV
+)
+volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten)
+left_boundary = F.SurfaceSubdomain1D(id=1, x=0)
+right_boundary = F.SurfaceSubdomain1D(id=2, x=L)
+
+my_model.subdomains = [
+    volume_subdomain,
+    left_boundary,
+    right_boundary,
+]
+```
+
+Now, we define our `incident_flux` and `gaussian_distribution` function. Here, we define the mean implantation depth $R_p$ and distribution width to be a few nanometers long:
+
+$$ R_p = 4 \times 10^{-9} \mathrm{m} $$
+$$ \mathrm{width} = 2.5 \times 10^{-9} \mathrm{m} $$
+
+```{code-cell} ipython3
+incident_flux = 1e13
+Rp = 4e-9
+width = 2.5e-9
+
+def gaussian_distribution(x, center, width):
+    return (
+        1
+        / (width * (2 * ufl.pi) ** 0.5)
+        * ufl.exp(-0.5 * ((x[0] - center) / width) ** 2)
+    )
+```
+
+We can define our species and use `ParticleSource` to represent the source term, and then add it to our model:
+
+```{code-cell} ipython3
+H = F.Species("H")
+my_model.species = [H]
+
+source_term = F.ParticleSource(
+    value=lambda x: incident_flux * gaussian_distribution(x, Rp, width),
+    volume=volume_subdomain,
+    species=H,
+)
+
+my_model.sources = [source_term]
+```
+
+Finally, we assign boundary conditions (zero concentration at $x=0$ and $x=L$) 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]
+
+plt.plot(x, c)
+plt.xlabel("x")
+plt.ylabel("c")
+plt.show()
+```
+
+We see that a huge spike in concentration in the first few nanometers of tungesten, where the implantation range is focused.
+
++++
+
+### Approximating plasma implantation using fixed concentration boundary conditions
+
+If recombination is fast enough, the spike shown above can be approximated as a fixed concentration boundary condition that mainly drives diffusion across the material. Learn more about the plasma implantation approximation approach _[here](https://festim.readthedocs.io/en/fenicsx/theory.html#plasma-implantation-approximation)_.
+
+To see how we might approximate this, let's define a maximum concentration to set on the left boundary, representing the spike from the implantation:
+
+$$ c_m = \frac{R_p \cdot \varphi_{imp}}{D} $$
+
+where $\varphi_{imp}$ is the implantation flux and $R_p$ is the implantation depth, both which we defined earlier, and $D$ is the material diffusivity:
+
+```{code-cell} ipython3
+D = tungsten.D_0 * np.exp(-tungsten.E_D / (F.k_B * my_model.temperature))
+c_m = incident_flux * Rp / D
+```
+
+Now, we'll change our boundary conditions to represent the implantation as a fixed concentration on the left boundary, and remove the source term from our problem:
+
+```{code-cell} ipython3
+my_model.boundary_conditions = [
+    F.FixedConcentrationBC(subdomain=left_boundary, value=c_m, species=H),
+    F.FixedConcentrationBC(subdomain=right_boundary, value=0, species=H),
+]
+
+my_model.sources = []
+profile_export = F.Profile1DExport(field=H,subdomain=volume_subdomain)
+my_model.exports = [profile_export]
+
+my_model.initialise()
+my_model.run()
+```
+
+Let's compare the profiles from the approximation to the volumetric source results:
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import matplotlib.pyplot as plt
+
+x2 = my_model.exports[0].x
+c2 = my_model.exports[0].data[0]
+
+plt.plot(x, c, label="With volumetric source")
+plt.plot(x2, c2, label="Approximation")
+
+plt.xlabel("x")
+plt.ylabel("c")
+plt.legend()
+plt.show()
+```
+
+Using the approximation is computationally less expensive, and still provides similar diffusion profiles.
+
++++
+
+## Complex isotopic exchange with multiple 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} \rightleftharpoons \mathrm{HD} \longrightarrow R_1 = K_{r1} c_H c_D - K_{d1}P_{HD} $$
+$$ \text{Reaction 2}: \mathrm{D+D} \rightleftharpoons \mathrm{D_2} \longrightarrow R_2 = K_{r2} c_D^2 - K_{d2}P_{D2} $$
+$$ \text{Reaction 3}: \mathrm{H+H} \rightleftharpoons \mathrm{H_2} \longrightarrow R_3 = K_{r3} c_H^2 - K_{d3}P_{H2} $$
+
+$$ \text{Reaction 4}: \mathrm{D+H_2} \rightleftharpoons \mathrm{HD + H} \longrightarrow R_4 = K_{r4} P_{H2} c_D - K_{d4}P_{HD} c_H $$
+
+
+From this reaction scheme, the surface fluxes of H and D are:
+
+$$
+\varphi_\mathrm{H} = -R_1  - 2 R_3 + R_4
+$$
+
+$$
+\varphi_\mathrm{D} = -R_1 - 2 R_2 - R_4
+$$
+
+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]
+
+P_h2 = 1000
+
+reaction_1_bc = 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,
+)
+
+reaction_2_bc = 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,
+)
+
+reaction_3_bc = F.SurfaceReactionBC(
+    reactant=[H, H],
+    gas_pressure=P_h2,
+    k_r0=1,
+    E_kr=0.1,
+    k_d0=1e-5,
+    E_kd=0.1,
+    subdomain=right_surf,
+)
+```
+
+Now, let's add our isotopic exchange reaction using `ParticleFluxBC` (see [here](h_transport_advanced.md) to learn more about defining isotopic exchange fluxes):
+
+```{code-cell} ipython3
+Kr_0_custom = 10000.0
+E_Kr_custom = 0.5  # eV
+
+def K_exchange(T):
+    return Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))
+
+def isotopic_exchange_D(c_D, T):
+    return - K_exchange(T) * P_h2 * c_D
+
+def isotopic_exchange_H(c_D, T):
+    return + K_exchange(T) * P_h2 * c_D
+
+
+reaction_4_D = F.ParticleFluxBC(
+    value=isotopic_exchange_D,
+    subdomain=right_surf,
+    species_dependent_value={"c_D": D},
+    species=D,
+)
+
+reaction_4_H = F.ParticleFluxBC(
+    value=isotopic_exchange_H,
+    subdomain=right_surf,
+    species_dependent_value={"c_D": D},
+    species=H,
+)
+```
+
+Finally, we add our boundary conditions and solve the steady-state problem:
+
+```{code-cell} ipython3
+my_model.boundary_conditions = [
+    reaction_1_bc,
+    reaction_2_bc,
+    reaction_3_bc,
+    reaction_4_D,
+    reaction_4_H,
+    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.
diff --git a/book/content/boundary_conditions/h_transport_advanced.md b/book/content/boundary_conditions/h_transport_advanced.md
new file mode 100644
index 0000000..7856375
--- /dev/null
+++ b/book/content/boundary_conditions/h_transport_advanced.md
@@ -0,0 +1,478 @@
+---
+jupytext:
+  formats: md:myst,ipynb
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.20.0
+kernelspec:
+  display_name: festim-workshop
+  language: python
+  name: python3
+---
+
+# Hydrogen transport: advanced
+
+This section discusses how to implement advanced boundary conditions (BCs) for hydrogen transport problems.
+
+Objectives:
+* Choose solubility laws (Henry's or Siervert's)
+* Implement surface reaction boundary conditions
+* Model isotopic exchange as a recombination flux
+
++++
+
+## Understanding background for solubility laws and surface reactions
+
++++
+
+### Solubility laws
+
+A solubility law prescribes the equilibrium relationship between the hydrogen concentration in a solid and the surrounding environment. Two common solubility laws for hydrogen are **Henry’s law** and **Sieverts’ law**:
+
+#### Henry's law
+
+Henry’s law assumes a linear relationship between the hydrogen concentration in the solid and the partial pressure of hydrogen in the gas:
+
+$$
+c_s = k_H \, P_H
+$$
+
+where:  
+- $c_s$ is the surface concentration,  
+- $P_H$ is the hydrogen partial pressure in the surrounding environment,  
+- $k_H$ is the Henry’s law constant (material and temperature-dependent).
+
+```{tip}
+Users interested in molten salt interactions with hydrogen will usually need to use *Henry's Law*.
+```
+
+#### Sieverts' law
+
+Sieverts’ law applies when hydrogen dissolves in metals as atomic hydrogen. The solubility is proportional to the square root of the hydrogen partial pressure (due to diatomic hydrogen dissociating into atoms in the solid):
+
+$$
+c_s = k_S \, \sqrt{P_H}
+$$
+
+where:  
+- $k_S$ is the Sieverts’ constant (which depends on temperature and the metal)
+
+```{tip}
+Users interested in liquid metal interactions with hydrogen will usually need to use **Sievert's Law**.
+```
+
++++
+
+### Imposing a solubility law
+
+Users can define the surface concentration using either Sieverts’ law. For Sieverts' law of solubility, we can use `SievertsBC`:
+
+```{code-cell} ipython3
+import festim as F
+
+boundary = F.SurfaceSubdomain(id=1)
+H = F.Species(name="Hydrogen")
+
+custom_pressure_value = lambda t: 2 + t
+my_sieverts_bc = F.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 `HenrysBC`:
+
+```{code-cell} ipython3
+pressure_value = lambda t: 5 * t
+my_henrys_bc = F.HenrysBC(subdomain=3, H_0=1.5, E_H=0.2, species=H, pressure=pressure_value)
+```
+
+To use these BCs in your simulation, add them to `boundary_conditions`:
+
+```{code-cell} ipython3
+my_model = F.HydrogenTransportProblem()
+my_model.boundary_conditions = [my_sieverts_bc, my_henrys_bc]
+```
+
+```{note}
+Using `HenrysBC` or `SievertsBC` is just a convenient way for users to define a `FixedConcentrationBC`, using the solubility law to calculate the concentration.
+```
+
++++
+
+---
+
++++
+
+## Simple surface reaction
+
+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}, \mathrm{B}$ are surface reactants  
+- $\mathrm{C}$ is the product species in the gas phase  
+- $\mathrm{K}_r$ and $\mathrm{K}_d$ are forward and backward reaction rates, respectively
+
+### Mathematical formulation
+
+Let:  
+- $c_A, c_B$ be the surface concentrations of $\mathrm{A}$ and $\mathrm{B}$  
+- $P_C$ be the partial pressure of $\mathrm{C}$  
+- $T$ the surface temperature  
+- $k_B$ the Boltzmann constant  
+
+The forward and backward reaction rates follow Arrhenius laws:
+
+$$
+K_r = k_{r0} \exp\!\left(-\frac{E_{kr}}{k_B T}\right), \qquad
+K_d = k_{d0} \exp\!\left(-\frac{E_{kd}}{k_B T}\right)
+$$
+
+The net surface reaction rate is:
+
+$$
+K = K_r \, c_A c_B - K_d \, P_C
+$$
+
+The resulting flux of species $\mathrm{A}$ into the surface is:
+
+$$
+\mathbf{J}_A \cdot \mathbf{n} = K
+$$
+
+If $\mathrm{A=B}$, the total flux becomes:
+
+$$
+\mathbf{J}_A \cdot \mathbf{n} = 2 K
+$$
+
+where $\mathbf{n}$ is the outward normal vector at the surface.
+
+### Recombination and dissociation
+
+If $\text{Species A = Species B}$, recombination or dissociation can also be modeled with `SurfaceReactionBC`, e.g.:
+
+$$
+\mathrm{H + H} \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad \mathrm{H_2}
+$$
+
+where 
+
+$$
+K = K_r \, c_H^2 - K_d \, P_{H_2}
+$$
+
+and corresponding flux of atomic hydrogen is:
+
+$$
+\mathbf{J}_H \cdot \mathbf{n} = - 2 K
+$$
+
++++
+
+### Implementing surface reaction boundary conditions
+
+Users can impose a surface reaction on a boundary using `SurfaceReactionBC`. 
+
+First, create a boundary to impose the BC and the reactant species:
+
+```{code-cell} ipython3
+import festim as F
+
+boundary = F.SurfaceSubdomain(id=1)
+A = F.Species("A")
+B = F.Species("B")
+```
+
+Now, we add these as arguments to the `SurfaceReactionBC` class. We'll also need to assign a `gas_pressure` (corresponding to the partial pressure of the product species), and the forward/backward rate (`k_r0`, `k_d0`) and energy (`E_kr`, `E_kd`) coefficients (see above to learn more about these parameters):
+
+```{code-cell} ipython3
+my_bc = F.SurfaceReactionBC(
+    reactant=[A, B],
+    gas_pressure=1e5,
+    k_r0=1,
+    E_kr=0.1,
+    k_d0=0,
+    E_kd=0,
+    subdomain=boundary,
+)
+```
+
+Finally, add the BC to your problem's `boundary_conditions` attribute using a list:
+
+```{code-cell} ipython3
+my_model = F.HydrogenTransportProblem()
+my_model.boundary_conditions = [my_bc]
+```
+
+```{note}
+Using `SurfaceReactionBC` is just a convenient way for users to define surface fluxes, which can be done manually (as shown below in the isotopic exchange example).
+```
+
++++
+
+---
+
++++
+
+## Isotopic exchange
+
+Isotopic exchange occurs when hydrogenic isotopes swap positions, for example:
+
+```{math}
+:label: t2_recomb
+\mathrm{T + T} \rightleftharpoons \mathrm{T_2}
+```
+
+```{math}
+:label: isotopic_exchange
+\mathrm{T + H_2} \rightleftharpoons \mathrm{H} + \mathrm{HT}
+```
+
+### Mathematical formulation
+
+Let:  
+- $c_T$ be the surface concentrations of tritium
+- $P_\mathrm{H_2}, P_\mathrm{HT}$ the partial pressures of $\mathrm{H_2}$ and $\mathrm{HT}$
+- $K_{r}$ the rate coefficient of reaction {eq}`t2_recomb`
+- $K_{r}^\ast$ the rate coefficient of reaction {eq}`isotopic_exchange`
+
+$$
+K_{r} = K_{r0} \exp\!\left(-\frac{E_{Kr}}{k_B T}\right)
+$$
+
+$$
+K_{r}^\ast = K_{r0}^\ast \exp\!\left(-\frac{E_{Kr}^\ast}{k_B T}\right)
+$$
+
+If $c_{H_2} \gg c_{HT}$, the flux of tritium is approximated as:
+
+$$
+\phi_T = \phi_\mathrm{recombination} + \phi_\mathrm{exchange}
+$$
+with
+$$
+\phi_\mathrm{recombination} = - K_{r} c_T^2 \\
+\phi_\mathrm{exchange} = - K_{r}^\ast P_{H_2} c_T
+$$
+
+These fluxes can be implemented in FESTIM using `ParticleFluxBC` with user-defined expressions for each reaction term, as shown below.
+
++++
+
+```{note}
+$K_{r0}^\ast$ and $K_{r0}$ have different dimensions since $P_\mathrm{H_2}$ is a pressure.
+```
+
+```{note}
+This example is not tracking the dissolution and transport of protium in the metal. And assumes that the partial presure of $\mathrm{H_2}$ is constant.
+```
+
++++
+
+```{note}
+For more complex isotopic exchanges, we can also use `SurfaceReactionBC` add other reactions. See [modeling isotopic exchange for multiple species](examples.md) to learn more.
+```
+
++++
+
+### Implementation
+
++++
+
+Let's work through a 1D example to illustrate the effect of isotopic exchange. We'll consider tritium diffusion from left to right, with a large partial pressure of $H_2$ at the right boundary (where isotopic exchange can occur).
+
+First, we'll set up our mesh and materials:
+
+```{code-cell} ipython3
+import numpy as np
+
+my_model = F.HydrogenTransportProblem()
+my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 100))
+
+right_subdomain = F.SurfaceSubdomain1D(id=2, x=1)
+left_subdomain = F.SurfaceSubdomain1D(id=3, x=0)
+
+material = F.Material(D_0=1e-2, E_D=0)
+
+vol = F.VolumeSubdomain1D(id=5, material=material, borders=[0, 1])
+my_model.subdomains = [vol, left_subdomain, right_subdomain]
+```
+
+Isotopic exchange reactions can be modeled as user-defined expressions using `ufl`. 
+
+```{code-cell} ipython3
+import ufl
+import festim as F
+
+Kr_0_custom = 1
+E_Kr_custom = 0.0  # eV
+P_H2 = 10  # assumed constant H2 concentration in
+
+def isotopic_exchange_flux_func(c, T):
+    return - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * P_H2 * c
+```
+
+Now we'll add our species `tritium` to our model and define a custom flux to represent isotopic exchange. We use `ParticleFluxBC` and the custom flux function we defined above to define this BC on the right boundary:
+
+```{code-cell} ipython3
+tritium = F.Species("T")
+my_model.species = [tritium]
+
+isotopic_exchange_flux = F.ParticleFluxBC(
+    value=isotopic_exchange_flux_func,
+    subdomain=right_subdomain,
+    species_dependent_value={"c": tritium},
+    species=tritium,
+)
+```
+
+```{code-cell} ipython3
+t2_recomb = F.SurfaceReactionBC(
+    reactant=[tritium, tritium],
+    gas_pressure=0,
+    k_r0=1e-22,
+    E_kr=0,
+    k_d0=0,
+    E_kd=0,
+    subdomain=right_subdomain,
+)
+```
+
+For diffusion across the mesh, we also define a fixed concentration on the left surface:
+
+```{code-cell} ipython3
+my_model.boundary_conditions = [
+    isotopic_exchange_flux,
+    t2_recomb,
+    F.FixedConcentrationBC(subdomain=left_subdomain, value=1e20, species=tritium),
+]
+```
+
+We'll export the flux using `SurfaceFlux` (see [exporting derived quantities to learn more](../post_process/derived.md)) and save this data to a variable `data_with_H2`. We also create a `Profile1DExport` to be able to plot the concentration profile later:
+
+```{code-cell} ipython3
+surface_flux = F.SurfaceFlux(field=tritium, surface=right_subdomain)
+
+profile = F.Profile1DExport(field=tritium, subdomain=vol)
+```
+
+```{code-cell} ipython3
+my_model.temperature = 300
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, stepsize=1, final_time=120)
+
+my_model.exports = [surface_flux, profile]
+my_model.initialise()
+my_model.run()
+
+data_with_H2 = surface_flux.data.copy()
+```
+
+Let's plot the evolution of the tritium concentration profile:
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import matplotlib.pyplot as plt
+
+import matplotlib.animation as animation
+from IPython.display import HTML
+
+# Turn off interactive plotting to prevent static display
+plt.ioff()
+
+# Create the figure and axis
+fig, ax = plt.subplots(figsize=(10, 6))
+
+# Set up the plot limits and labels
+ax.set_ylabel("Tritium concentration")
+ax.set_xlabel("Position")
+
+# Initialize empty line objects for the animation
+line_current, = ax.plot([], [], color="black", linewidth=2, zorder=1000)
+lines_previous = []
+
+def animate(frame):
+    # Clear previous iteration lines
+    for line in lines_previous:
+        line.remove()
+    lines_previous.clear()
+    
+    # Plot all previous iterations in light grey
+    for i in range(frame):
+        line_prev, = ax.plot(profile.x, profile.data[i], color="lightgrey", linewidth=0.8, alpha=0.6)
+        lines_previous.append(line_prev)
+    
+    # Plot current iteration in black
+    line_current.set_data(profile.x, profile.data[frame])
+    
+    ax.set_title(f"Time: {profile.t[frame]:.2f}")
+
+    return [line_current] + lines_previous
+
+# Create animation
+anim = animation.FuncAnimation(
+    fig, animate, frames=len(profile.data), interval=100, blit=False, repeat=True
+)
+
+# Close the figure to prevent static display
+plt.close(fig)
+
+# Display only the animation
+HTML(anim.to_jshtml())
+```
+
+We can compare the surface flux to the case where we have no isotopic exchange by removing the isotopic exchange boundary condition and saving the results to `data_no_H2`:
+
+```{code-cell} ipython3
+my_model.boundary_conditions = [
+    F.FixedConcentrationBC(subdomain=left_subdomain, value=1e20, species=tritium),
+    t2_recomb,
+]
+
+profile = F.Profile1DExport(field=tritium, subdomain=vol)
+my_model.exports = [surface_flux, profile]
+my_model.initialise()
+my_model.run()
+data_no_H2 = surface_flux.data.copy()
+```
+
+Without isotopic exchange the gradient is much lower, indicating that isotopic exchange helps enhance surface outgassing of tritium!
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+# Create animation
+anim = animation.FuncAnimation(
+    fig, animate, frames=len(profile.data), interval=100, blit=False, repeat=True
+)
+
+# Close the figure to prevent static display
+plt.close(fig)
+
+# Display only the animation
+HTML(anim.to_jshtml())
+```
+
+Let's plot both `data_with_H2` and `data_no_H2` to see the temporal evolution of the surface flux for each case:
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import matplotlib.pyplot as plt
+
+x = my_model.mesh.mesh.geometry.x[:, 0]
+t = surface_flux.t
+plt.plot(t, data_no_H2, label="No H_2")
+plt.plot(t, data_with_H2,label="With H_2")
+
+plt.xlabel("Time")
+plt.ylabel("Surface Flux (right)")
+plt.legend(reverse=True)
+plt.show()
+```
+
+We see that the flux is higher with isotopic exchange, as we'd expect.
diff --git a/book/content/boundary_conditions/h_transport_basic.md b/book/content/boundary_conditions/h_transport_basic.md
new file mode 100644
index 0000000..4ca7d6a
--- /dev/null
+++ b/book/content/boundary_conditions/h_transport_basic.md
@@ -0,0 +1,401 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.19.1
+kernelspec:
+  display_name: festim-workshop
+  language: python
+  name: python3
+---
+
+# Hydrogen transport: basic
+
+This section discusses how to implement basic boundary conditions (BCs) for hydrogen transport problems. Boundary conditions 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)_.
+
+Objectives:
+* Implement fixed concentration boundary conditions
+* Implement particle flux boundary conditions
+
++++
+
+## Imposing fixed concentration boundary conditions
+
+Users can prescribe a fixed concentration on a boundary using `FixedConcentrationBC`, which can depend on temperature, time, and space.
+
+Let us consider a 2D, steady-state, single-species example with the following *space-dependent* boundary conditions:
+
+$$ \text{Top:} \quad c = 2x + y $$
+$$ \text{Right:} \quad c = y^2 + 1 $$
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import festim as F
+import numpy as np
+from dolfinx.mesh import create_unit_square
+from mpi4py import MPI
+
+my_model = F.HydrogenTransportProblem()
+my_model.mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, 10, 10))
+H = F.Species("H")
+my_model.species = [H]
+mat = F.Material(D_0=1e-3, E_D=0)
+my_model.material = mat
+my_model.temperature = 400
+my_model.settings = F.Settings(atol=1e-8, rtol=1e-8, transient=False)
+```
+
+Now, let's create the boundaries to assign the BCs:
+
+```{code-cell} ipython3
+top_subdomain = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1))
+right_subdomain = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1))
+left_subdomain = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0))
+bottom_subdomain = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[1], 0))
+vol = F.VolumeSubdomain(id=5, material=mat)
+
+my_model.subdomains = [top_subdomain, right_subdomain, left_subdomain, bottom_subdomain, vol]
+```
+
+We can define the top and right boundary conditions using `lambda` functions:
+
+```{code-cell} ipython3
+top_bc_function = lambda x: 2.0 * x[0] + x[1]
+right_bc_function = lambda x: x[1] ** 2 + 1.0
+```
+
+```{note}
+`x[0]`, `x[1]`, `x[2]` correspond to the first, second, and third space coordinate ( *(x, y, z)* in cartesian space). See *{ref}`label-time-temp-concentration`* to learn more.
+```
+
++++
+
+To include these boundary conditions to our problem, we use `FixedConcentrationBC`. We must also specify which subdomain (boundary) each BC is applied to, as well as the corresponding species:
+
+```{code-cell} ipython3
+top_bc = F.FixedConcentrationBC(
+    subdomain=top_subdomain,
+    value=top_bc_function,
+    species=H
+)
+
+right_bc = F.FixedConcentrationBC(
+    subdomain=right_subdomain,
+    value=right_bc_function,
+    species=H
+)
+```
+
+We can add these to our problem, we pass them into `boundary_conditions` as a list and then run the simulation:
+
+```{code-cell} ipython3
+my_model.boundary_conditions = [top_bc, right_bc]
+my_model.initialise()
+my_model.run()
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import pyvista
+from dolfinx import plot
+
+topology, cell_types, geometry = plot.vtk_mesh(H.post_processing_solution.function_space)
+grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
+grid.point_data["c"] = H.post_processing_solution.x.array
+grid.set_active_scalars("c")
+
+pyvista.set_jupyter_backend("html")
+
+plotter = pyvista.Plotter()
+
+plotter.add_mesh(grid)
+plotter.view_xy()
+contours = grid.contour(isosurfaces=5, scalars="c")
+if not pyvista.OFF_SCREEN:
+    plotter.show()
+else:
+    figure = plotter.screenshot("concentration.png")
+```
+
+```{note}
+
+In this example, we did not define any boundary conditions for the left and bottom surfaces. If no BC is set on a boundary, it is assumed that the flux is null (symmetry boundary condition):
+
+$$
+\mathbf{J} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma
+$$
+
+```
+
++++
+
+(label-time-temp-concentration)=
+### Adding time and temperature dependent boundary conditions ###
+
++++
+
+Users can also impose concentration BCs that are dependent on space $x$, time $t$, and temperature $T$, such as:
+
+$$ c = 10 + x^2 + T(t+2) $$
+
+To do so, we define a Python function:
+
+```{code-cell} ipython3
+def my_custom_value(x, t, T):
+    return 10 + x[0]**2 + T *(t + 2)
+
+boundary = F.SurfaceSubdomain(id=1)
+my_bc = F.FixedConcentrationBC(subdomain=boundary, value=my_custom_value, species=H)
+```
+
+```{tip}
+This custom function is equivalent to:
+
+`
+my_custom_value = lambda x, t, T: 10 + x[0]**2 + T *(t + 2)
+`
+```
+
+### Multi-species fixed concentration BC ##
+
+Users may also want to add BCs for multiple species in their model. 
+
+For example, if we wanted to impose separate fixed concentrations for deuterium and hydrogen on a boundary, we need to define each species and then one BC for each:
+
+```{code-cell} ipython3
+import festim as F
+
+H = F.Species("H")
+D = F.Species("D")
+boundary = F.SurfaceSubdomain(id=1)
+
+top_H = F.FixedConcentrationBC(
+    subdomain=boundary,
+    value=5.0,
+    species=H
+)
+
+top_D = F.FixedConcentrationBC(
+    subdomain=boundary,
+    value=10.0,
+    species=D
+)       
+```
+
+## Imposing particle flux boundary conditions
+
++++
+
+Users can impose the particle flux at boundaries using `ParticleFluxBC` class. At minimum, this BC requires a boundary, value, and species to be added:
+
+```{code-cell} ipython3
+import festim as F
+
+boundary = F.SurfaceSubdomain(id=1)
+H = F.Species(name="H")
+
+my_flux_bc = F.ParticleFluxBC(subdomain=boundary, value=2, species=H)
+```
+
+To use this BC to your problem, simply add it to `boundary_conditions` attribute of your problem:
+
+```{code-cell} ipython3
+my_model = F.HydrogenTransportProblem()
+my_model.boundary_conditions = [my_flux_bc]
+```
+
+### Species-dependent 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
+import festim as F
+
+my_model = F.HydrogenTransportProblem()
+boundary = F.SurfaceSubdomain(id=1)
+H = F.Species(name="H")
+
+J = lambda t, c: 10*t**2 + 2*c
+```
+
+```{note}
+A flux boundary condition dependent on the actual solution is called a Robin boundary condition.
+```
+
++++
+
+To add this BC, we need to create a dictionary that maps the concentration variable `c` in our custom function to our species `H`:
+
+```{code-cell} ipython3
+mapping_dict = {"c": H}
+```
+
+Now, we add our `ParticleFluxBC`, also utilizing the `species_dependent_value` argument:
+
+```{code-cell} ipython3
+my_flux_bc = F.ParticleFluxBC(
+    subdomain=boundary,
+    value=J,
+    species=H,
+    species_dependent_value=mapping_dict,
+)
+
+my_model.boundary_conditions = [my_flux_bc]
+```
+
+### 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}$. The boundary conditions include recombination fluxes $\phi_{AB}$ and $\phi_{ABC}$ that depend on the concentrations $\mathrm{c}$ (on the right) and fixed concentrations (on the left):
+
+$$ \text{Right (recombination flux):} \quad \phi_{AB} = -c_A c_B $$
+
+$$ \text{Right (recombination flux):} \quad \phi_{ABC} = -10 c_A c_B c_C$$
+
+$$ \text{Left (fixed concentration):} \quad c_A = c_B = c_C = 1 $$
+
+
+First, let us define each species using `Species`:
+
+```{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]
+```
+
+Now, we create the dictionary to be passed into `species_dependent_value`; this maps each argument in the custom flux function to the corresponding defined species:
+
+```{code-cell} ipython3
+species_dependent_value = {"c_A": A, "c_B": B, "c_C": C}
+```
+
+We also define our custom functions for the fluxes:
+
+```{code-cell} ipython3
+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 corresponding boundaries:
+
+```{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]
+```
+
+```{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](../material/material_advanced.md). Here we choose arbitrarily different diffusivity coefficients to see transport between species. 
+```
+
++++
+
+Let's 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
+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),     
+]
+```
+
+We can export the flux for each species on the right using `SurfaceFlux` (see [post-processing derived values](../post_process/derived.md) to learn more about exporting fluxes):
+
+```{code-cell} ipython3
+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,
+]
+```
+
+Finally, let's solve and plot the profile for each species:
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+my_model.temperature = 300
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
+
+my_model.initialise()
+my_model.run()
+
+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/environment.yml b/environment.yml
index e021c3f..bed85ab 100644
--- a/environment.yml
+++ b/environment.yml
@@ -2,7 +2,12 @@ name: festim-workshop
 channels:
   - conda-forge
   - defaults
+
 dependencies:
+  - python<3.14
+  - pip=23.3
+  - setuptools<82
+  - wheel<0.45
   - jupyter-book<2.0
   - python<3.14  # som pyvista incompatibility issues with python 3.14
   - pip=23.3