diff --git a/test/tests/val-flibe/run_val_flibe.py b/test/tests/val-flibe/run_val_flibe.py
new file mode 100644
index 00000000..1211f780
--- /dev/null
+++ b/test/tests/val-flibe/run_val_flibe.py
@@ -0,0 +1,48 @@
+import numpy as np
+import subprocess
+
+# number of processors
+n_proc = 4
+
+filename = 'val-flibe.i'
+search_strings = ["T = '${units ",'p_bnd = ','    file_base = ']
+command_to_run = "mpirun -np "+str(n_proc)+" ~/projects/TMAP8/tmap8-opt -i val-flibe.i"
+
+# Experimental data pressure and Temperature series
+p_exp = np.array([1210,538,315,171,1210,538,315,1210,538,1210])
+T_exp = np.array([700,700,700,700,650,650,650,600,600,550])
+
+
+def run_terminal_command(command):
+    subprocess.run([command], shell=True, check=True)
+
+def substitute_row(filename, search_string, replacement):
+    # Read the contents of the file
+    with open(filename, 'r') as file:
+        lines = file.readlines()
+    # Search for the row and substitute it
+    for i, line in enumerate(lines):
+        if search_string in line:
+            lines[i] = replacement + '\n'  # Add a newline character for consistency
+    # Write the modified contents back to the file
+    with open(filename, 'w') as file:
+        file.writelines(lines)
+
+
+
+
+# Names of output files for each case
+output_fname = np.array([])
+for p,T in zip(p_exp,T_exp):
+    output_fname = np.append(output_fname, 'val-flibe_'+str(p)+'_'+str(T))
+
+
+
+for T,p,outName in zip(T_exp,p_exp,output_fname):
+    replacemens = ["T = '${units "+str(T)+" degC -> K}' # temperature", 'p_bnd = '+str(p)+' # pressure', "    file_base = '"+outName+"'"]
+    substitute_row(filename, search_strings[0], replacemens[0])
+    substitute_row(filename, search_strings[1], replacemens[1])
+    substitute_row(filename, search_strings[2], replacemens[2])
+    # Run input file with new BCs
+    run_terminal_command(command_to_run)
+    print('########## Run at',p,'Pa and',T,'degC completed. ##########')
diff --git a/test/tests/val-flibe/tests b/test/tests/val-flibe/tests
new file mode 100644
index 00000000..856b83e3
--- /dev/null
+++ b/test/tests/val-flibe/tests
@@ -0,0 +1,11 @@
+[Tests]
+  design = 'ADMatDiffusion.md TimeDerivative.md InterfaceSorption.md EquilibriumBC.md ADDirichletBC.md'
+  issues = '#193'
+    [val-flibe]
+      type = CSVDiff
+      input = val-flibe.i
+	    should_execute = True
+	    csvdiff = val-flibe_out.csv
+      requirement = 'The system shall be able to model transient diffusion through a crucible containing molten salt with a constant concentration boundary condition as the species source.'
+    []
+[]
diff --git a/test/tests/val-flibe/val-flibe.i b/test/tests/val-flibe/val-flibe.i
new file mode 100644
index 00000000..1f47ee4b
--- /dev/null
+++ b/test/tests/val-flibe/val-flibe.i
@@ -0,0 +1,285 @@
+endtime = 140000 # simulation end time
+
+R = 8.31446261815324 # Gas constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)
+
+T = '${units 550 degC -> K}' # temperature
+
+p_bnd = 1210 # pressure
+
+L_Ni = '${units 2 mm -> m}' # nickel thickness
+L_salt = '${units 8.1 mm -> m}' # salt thickness
+
+num_nodes = 200 # (-)
+
+[Mesh]
+  [whole_domain]
+    type = GeneratedMeshGenerator
+    xmin = 0
+    xmax = '${fparse L_Ni + L_salt}'
+    dim = 1
+    nx = ${num_nodes}
+  []
+  [block_1]
+    type = ParsedSubdomainMeshGenerator
+    input = whole_domain
+    combinatorial_geometry = 'x <= ${L_Ni}'
+    block_id = 0
+  []
+  [block_2]
+    type = ParsedSubdomainMeshGenerator
+    input = block_1
+    combinatorial_geometry = 'x > ${L_Ni}'
+    block_id = 1
+  []
+  [interface]
+    type = SideSetsBetweenSubdomainsGenerator
+    input = block_2
+    primary_block = '0' # Ni
+    paired_block = '1' # salt
+    new_boundary = 'interface'
+  []
+  [interface_other_side]
+    type = SideSetsBetweenSubdomainsGenerator
+    input = interface
+    primary_block = '1' # salt
+    paired_block = '0' # Ni
+    new_boundary = 'interface_other'
+  []
+[]
+
+[Variables]
+  [conc_Ni]
+    initial_condition = 0.9e-0
+    block = 0
+  []
+  [conc_salt]
+    initial_condition = 0.3e-0
+    block = 1
+  []
+[]
+
+[AuxVariables]
+  [enclosure_pressure]
+    family = SCALAR
+    initial_condition = ${p_bnd}
+  []
+  [flux_x]
+    order = FIRST
+    family = MONOMIAL
+  []
+  [p_div_RT_salt]
+    block = 1
+    initial_condition = 0.001
+  []
+  [p_div_RT_Ni]
+    block = 0
+    initial_condition = 0.001
+  []
+[]
+
+[Kernels]
+  [diff_Ni]
+    type = ADMatDiffusion
+    variable = conc_Ni
+    diffusivity = diffusivity_Ni
+    block = 0
+  []
+  [diff_salt]
+    type = ADMatDiffusion
+    variable = conc_salt
+    diffusivity = diffusivity_salt
+    block = 1
+  []
+  [time_diff_Ni]
+    type = TimeDerivative
+    variable = conc_Ni
+    block = 0
+  []
+  [time_diff_salt]
+    type = TimeDerivative
+    variable = conc_salt
+    block = 1
+  []
+[]
+
+[InterfaceKernels]
+  [tied]
+    type = InterfaceSorption
+    K0 = 0.564
+    Ea = 15800.0
+    n_sorption = 0.5
+    diffusivity = diffusivity_Ni_nonAD
+    unit_scale = 1
+    unit_scale_neighbor = 1
+    temperature = ${T}
+    variable = conc_Ni
+    neighbor_var = p_div_RT_salt
+    sorption_penalty = 0.1
+    boundary = 'interface'
+  []
+[]
+
+[AuxKernels]
+  [flux_x_Ni]
+    type = DiffusionFluxAux
+    diffusivity = diffusivity_Ni
+    variable = flux_x
+    diffusion_variable = conc_Ni
+    component = x
+    block = 0
+  []
+  [flux_x_salt]
+    type = DiffusionFluxAux
+    diffusivity = diffusivity_salt
+    variable = flux_x
+    diffusion_variable = conc_salt
+    component = x
+    block = 1
+  []
+  [p_div_RT_Ni_kernel]
+    variable = p_div_RT_Ni
+    type = ParsedAux
+    expression = '(conc_Ni / (${R} * ${T} * 0.564 * exp(-15800/(${R}*${T}))))^2'
+    coupled_variables = 'conc_Ni'
+  []
+  [p_div_RT_salt_kernel]
+    variable = p_div_RT_salt
+    type = ParsedAux
+    expression = '(conc_salt / (${R} * ${T} * 0.079 * exp(-35000/(${R}*${T}))))'
+    coupled_variables = 'conc_salt'
+  []
+[]
+
+[BCs]
+  [left_flux]
+    type = EquilibriumBC
+    Ko = 0.564
+    activation_energy = 15800.0
+    boundary = left
+    enclosure_var = enclosure_pressure
+    temperature = ${T}
+    variable = conc_Ni
+    p = 0.5
+  []
+  [right_flux]
+    type = ADDirichletBC
+    boundary = right
+    variable = conc_salt
+    value = 0.0
+  []
+[]
+
+[Functions]
+  [diffusivity_Ni_func]
+    type = ParsedFunction
+    symbol_names = 'T'
+    symbol_values = '${T}'
+    expression = '0.0000007*exp(-39500/(${R}*T))'
+  []
+
+  [diffusivity_salt_func]
+    type = ParsedFunction
+    symbol_names = 'T'
+    symbol_values = '${T}'
+    expression = '0.00000093*exp(-42000/(${R}*T))'
+  []
+
+  [solubility_Ni_func]
+    type = ParsedFunction
+    symbol_names = 'T'
+    symbol_values = '${T}'
+    expression = '0.564 * exp(-15800/(${R}*T))'
+  []
+
+  [solubility_salt_func]
+    type = ParsedFunction
+    symbol_names = 'T'
+    symbol_values = '${T}'
+    expression = '0.079 * exp(-35000/(${R}*T))'
+  []
+[]
+
+[Materials]
+  [diff_solu]
+    type = ADGenericFunctionMaterial
+    prop_names = 'diffusivity_Ni diffusivity_salt solubility_Ni solubility_salt'
+    prop_values = 'diffusivity_Ni_func diffusivity_salt_func solubility_Ni_func solubility_salt_func'
+    outputs = all
+  []
+  [converter_to_regular]
+    type = MaterialADConverter
+    ad_props_in = 'diffusivity_Ni diffusivity_salt'
+    reg_props_out = 'diffusivity_Ni_nonAD diffusivity_salt_nonAD'
+  []
+[]
+
+[Postprocessors]
+  [avg_flux_right]
+    type = SideDiffusiveFluxAverage
+    variable = conc_salt
+    boundary = right
+    diffusivity = diffusivity_salt_nonAD
+  []
+[]
+
+[Preconditioning]
+  [SMP]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  steady_state_detection = true
+  steady_state_start_time = 40000
+  steady_state_tolerance = 1e-9
+  scheme = bdf2 # bdf2 # crank-nicolson # explicit-euler
+  solve_type = NEWTON # LINEAR # JFNK # NEWTON
+  petsc_options = '-snes_ksp_ew'
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+  l_max_its = 10
+  nl_max_its = 13
+  nl_rel_tol = 1e-8 # nonlinear relative tolerance
+  nl_abs_tol = 1e-20 #1e-30 # nonlinear absolute tolerance
+  l_tol = 1e-6 # 1e-3 - 1e-5 # linear tolerance
+  end_time = ${endtime}
+  automatic_scaling = true
+  line_search = none
+  dtmax = 10.0
+  [TimeStepper]
+    type = IterationAdaptiveDT
+    dt = 1e-10
+    optimal_iterations = 18 # 6-10 or 18
+    growth_factor = 1.1
+    cutback_factor = 0.5
+  []
+[]
+
+[Outputs]
+  execute_on = timestep_end
+  exodus = true
+  [csv]
+    type = CSV
+    file_base = 'val-flibe_1210_550'
+    time_step_interval = 500
+  []
+[]
+
+[Dampers]
+  [limit_salt]
+    type = BoundingValueElementDamper
+    variable = conc_salt
+    max_value = 1e42
+    min_value = -0.01
+    min_damping = 0.001
+  []
+  [limit_Ni]
+    type = BoundingValueElementDamper
+    variable = conc_Ni
+    max_value = 1e42
+    min_value = -0.01
+    min_damping = 0.001
+  []
+[]
diff --git a/test/tests/val-flibe/val-flibe_test_BC.i b/test/tests/val-flibe/val-flibe_test_BC.i
new file mode 100644
index 00000000..b1b73acf
--- /dev/null
+++ b/test/tests/val-flibe/val-flibe_test_BC.i
@@ -0,0 +1,288 @@
+endtime = 140000 # simulation end time
+
+R = 8.31446261815324 # Gas constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)
+
+T = '${units 550 degC -> K}' # temperature
+
+p_bnd = 1210 # pressure
+
+L_Ni = '${units 2 mm -> m}' # nickel thickness
+L_salt = '${units 8.1 mm -> m}' # salt thickness
+
+num_nodes = 200 # (-)
+
+[Mesh]
+  [whole_domain]
+    type = GeneratedMeshGenerator
+    xmin = 0
+    xmax = '${fparse L_Ni + L_salt}'
+    dim = 1
+    nx = ${num_nodes}
+  []
+  [block_1]
+    type = ParsedSubdomainMeshGenerator
+    input = whole_domain
+    combinatorial_geometry = 'x <= ${L_Ni}'
+    block_id = 0
+  []
+  [block_2]
+    type = ParsedSubdomainMeshGenerator
+    input = block_1
+    combinatorial_geometry = 'x > ${L_Ni}'
+    block_id = 1
+  []
+  [interface]
+    type = SideSetsBetweenSubdomainsGenerator
+    input = block_2
+    primary_block = '0' # Ni
+    paired_block = '1' # salt
+    new_boundary = 'interface'
+  []
+  [interface_other_side]
+    type = SideSetsBetweenSubdomainsGenerator
+    input = interface
+    primary_block = '1' # salt
+    paired_block = '0' # Ni
+    new_boundary = 'interface_other'
+  []
+[]
+
+[Variables]
+  [conc_Ni]
+    initial_condition = 0.9e-0
+    block = 0
+  []
+  [conc_salt]
+    initial_condition = 0.3e-0
+    block = 1
+  []
+[]
+
+[AuxVariables]
+  [enclosure_pressure]
+    family = SCALAR
+    initial_condition = ${p_bnd}
+  []
+  [flux_x]
+    order = FIRST
+    family = MONOMIAL
+  []
+  [p_div_RT_salt]
+    block = 1
+    initial_condition = 0.001
+  []
+  [p_div_RT_Ni]
+    block = 0
+    initial_condition = 0.001
+  []
+[]
+
+[Kernels]
+  [diff_Ni]
+    type = ADMatDiffusion
+    variable = conc_Ni
+    diffusivity = diffusivity_Ni
+    block = 0
+  []
+  [diff_salt]
+    type = ADMatDiffusion
+    variable = conc_salt
+    diffusivity = diffusivity_salt
+    block = 1
+  []
+  [time_diff_Ni]
+    type = TimeDerivative
+    variable = conc_Ni
+    block = 0
+  []
+  [time_diff_salt]
+    type = TimeDerivative
+    variable = conc_salt
+    block = 1
+  []
+[]
+
+[InterfaceKernels]
+  [tied]
+    type = ADPenaltyInterfaceDiffusion
+    variable = conc_Ni
+    neighbor_var = conc_salt
+    penalty = 0.06 # to avoid cutting time-step on civet testing
+    jump_prop_name = solubility_ratio
+    boundary = 'interface'
+  []
+[]
+
+[AuxKernels]
+  [flux_x_Ni]
+    type = DiffusionFluxAux
+    diffusivity = diffusivity_Ni
+    variable = flux_x
+    diffusion_variable = conc_Ni
+    component = x
+    block = 0
+  []
+  [flux_x_salt]
+    type = DiffusionFluxAux
+    diffusivity = diffusivity_salt
+    variable = flux_x
+    diffusion_variable = conc_salt
+    component = x
+    block = 1
+  []
+  [p_div_RT_Ni_kernel]
+    variable = p_div_RT_Ni
+    type = ParsedAux
+    expression = '(conc_Ni / (${R} * ${T} * 0.564 * exp(-15800/(${R}*${T}))))^2'
+    coupled_variables = 'conc_Ni'
+  []
+  [p_div_RT_salt_kernel]
+    variable = p_div_RT_salt
+    type = ParsedAux
+    expression = '(conc_salt / (${R} * ${T} * 0.079 * exp(-35000/(${R}*${T}))))'
+    coupled_variables = 'conc_salt'
+  []
+[]
+
+[BCs]
+  [left_flux]
+    type = EquilibriumBC
+    Ko = 0.564
+    activation_energy = 15800.0
+    boundary = left
+    enclosure_var = enclosure_pressure
+    temperature = ${T}
+    variable = conc_Ni
+    p = 0.5
+  []
+  [right_flux]
+    type = ADDirichletBC
+    boundary = right
+    variable = conc_salt
+    value = 0.0
+  []
+[]
+
+[Functions]
+  [diffusivity_Ni_func]
+    type = ParsedFunction
+    symbol_names = 'T'
+    symbol_values = '${T}'
+    expression = '0.0000007*exp(-39500/(${R}*T))'
+  []
+
+  [diffusivity_salt_func]
+    type = ParsedFunction
+    symbol_names = 'T'
+    symbol_values = '${T}'
+    expression = '0.00000093*exp(-42000/(${R}*T))'
+  []
+
+  [solubility_Ni_func]
+    type = ParsedFunction
+    symbol_names = 'T'
+    symbol_values = '${T}'
+    expression = '0.564 * exp(-15800/(${R}*T))'
+  []
+
+  [solubility_salt_func]
+    type = ParsedFunction
+    symbol_names = 'T'
+    symbol_values = '${T}'
+    expression = '0.079 * exp(-35000/(${R}*T))'
+  []
+[]
+
+[Materials]
+  [diff_solu]
+    type = ADGenericFunctionMaterial
+    prop_names = 'diffusivity_Ni diffusivity_salt solubility_Ni solubility_salt'
+    prop_values = 'diffusivity_Ni_func diffusivity_salt_func solubility_Ni_func solubility_salt_func'
+    outputs = all
+  []
+  [converter_to_regular]
+    type = MaterialADConverter
+    ad_props_in = 'diffusivity_Ni diffusivity_salt'
+    reg_props_out = 'diffusivity_Ni_nonAD diffusivity_salt_nonAD'
+  []
+
+  [interface_jump]
+    type = SolubilityRatioMaterial
+    solubility_primary = solubility_Ni
+    solubility_secondary = solubility_salt
+    boundary = interface
+    concentration_primary = conc_Ni
+    concentration_secondary = conc_salt
+  []
+[]
+
+[Postprocessors]
+  [avg_flux_right]
+    type = SideDiffusiveFluxAverage
+    variable = conc_salt
+    boundary = right
+    diffusivity = diffusivity_salt_nonAD
+  []
+[]
+
+[Preconditioning]
+  [SMP]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  steady_state_detection = true
+  steady_state_start_time = 40000
+  steady_state_tolerance = 1e-9
+  scheme = bdf2 # bdf2 # crank-nicolson # explicit-euler
+  solve_type = NEWTON # LINEAR # JFNK # NEWTON
+  petsc_options = '-snes_ksp_ew'
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+  l_max_its = 10
+  nl_max_its = 13
+  nl_rel_tol = 1e-8 # nonlinear relative tolerance
+  nl_abs_tol = 1e-20 #1e-30 # nonlinear absolute tolerance
+  l_tol = 1e-6 # 1e-3 - 1e-5 # linear tolerance
+  end_time = ${endtime}
+  automatic_scaling = true
+  line_search = none
+  dtmax = 10.0
+  [TimeStepper]
+    type = IterationAdaptiveDT
+    dt = 1e-10
+    optimal_iterations = 18 # 6-10 or 18
+    growth_factor = 1.1
+    cutback_factor = 0.5
+  []
+[]
+
+[Outputs]
+  execute_on = timestep_end
+  exodus = true
+  [csv]
+    type = CSV
+    file_base = 'val-flibe_1210_550'
+    time_step_interval = 500
+  []
+[]
+
+[Dampers]
+  [limit_salt]
+    type = BoundingValueElementDamper
+    variable = conc_salt
+    max_value = 1e42
+    min_value = -0.01
+    min_damping = 0.001
+  []
+  [limit_Ni]
+    type = BoundingValueElementDamper
+    variable = conc_Ni
+    max_value = 1e42
+    min_value = -0.01
+    min_damping = 0.001
+  []
+[]