diff --git a/FabNESO/__init__.py b/FabNESO/__init__.py
index c5bfed2..9e82704 100644
--- a/FabNESO/__init__.py
+++ b/FabNESO/__init__.py
@@ -1,9 +1,23 @@
 """Neptune Exploratory SOftware (NESO) plugin for FabSim3."""
 
 try:
-    from .tasks import neso, neso_ensemble, neso_vbmc, neso_write_field
+    from .tasks import (
+        neso,
+        neso_ensemble,
+        neso_vbmc,
+        neso_vvuq_campaign,
+        neso_write_field,
+        single_run_vvuq,
+    )
 
-    __all__ = ["neso", "neso_ensemble", "neso_vbmc", "neso_write_field"]
+    __all__ = [
+        "neso",
+        "neso_ensemble",
+        "neso_vbmc",
+        "neso_write_field",
+        "neso_vvuq_campaign",
+        "single_run_vvuq",
+    ]
 
 except ImportError:
     import warnings
diff --git a/FabNESO/run_vvuq_instance.sh b/FabNESO/run_vvuq_instance.sh
new file mode 100755
index 0000000..5209d6e
--- /dev/null
+++ b/FabNESO/run_vvuq_instance.sh
@@ -0,0 +1,5 @@
+#!/usr/bin/env bash
+
+echo fabsim $1 single_run_vvuq:$PWD,mesh_file=$3,solver=$4,processes=$5,nodes=$6,cpus_per_process=$7,wall_time=$8,neso_outfile=$9
+
+fabsim $1 single_run_vvuq:$PWD,mesh_file=$3,solver=$4,processes=$5,nodes=$6,cpus_per_process=$7,wall_time=$8,neso_outfile=$9
diff --git a/FabNESO/tasks.py b/FabNESO/tasks.py
index 611441a..b023c30 100644
--- a/FabNESO/tasks.py
+++ b/FabNESO/tasks.py
@@ -7,11 +7,14 @@
 import re
 import shutil
 import time
+from ast import literal_eval
 from contextlib import nullcontext
 from pathlib import Path
 from tempfile import TemporaryDirectory
 from typing import Any
 
+import chaospy as cp
+import easyvvuq as uq
 import numpy as np
 import pyvbmc
 
@@ -579,3 +582,199 @@ def log_density(
         (config_dict["initial_run"] - observed_results["field_value"]) ** 2
         / (2 * config_dict["observation_noise_std"] ** 2)
     ).sum()
+
+
+def _parse_vvuq_parameters(
+    vvuq_param_file: Path,
+) -> tuple[dict[Any, dict[str, float]], dict[str, Any], list[str]]:
+    # read the inpurt file
+    with Path.open(vvuq_param_file) as f:
+        data = f.read()
+    d = literal_eval(data)
+    vvuq_params = {}
+    parameter_variations = {}
+    for parameter, param_values in d["parameters"].items():
+        vvuq_params[parameter] = {
+            "type": "float",
+            "default": param_values[0],
+        }
+        vary_length = 3
+        if len(param_values) == vary_length:
+            # Make the chaospy distribution. Possibly to be able to edit this later?
+            parameter_variations[parameter] = cp.Uniform(
+                param_values[1], param_values[2]
+            )
+    output_columns = d["output_columns"]
+    return (vvuq_params, parameter_variations, output_columns)
+
+
+@fab.task
+@fab.load_plugin_env_vars("FabNESO")
+def single_run_vvuq(
+    config: str,
+    mesh_file: str = "mesh.xml",
+    solver: str = "Electrostatic2D3V",
+    processes: str | int = 4,
+    nodes: str | int = 1,
+    cpus_per_process: str | int = 1,
+    wall_time: str = "00:15:00",
+    neso_outfile: str = "Electrostatic2D3V_line_field_evaluations.h5part",
+) -> None:
+    """
+    Run a single instance of the VVUQ campaign.
+
+    Run an instance of the VVUQ campaign using configurations in the designated solver.
+    The conditions and mesh are, by design, the defaults. When the iteration is done,
+    copy the NESO output to the local directory for later decoding through VVUQ.
+
+
+    Args:
+        config: Directory with single run configuration information.
+
+    Keyword Args:
+        solver: Which NESO solver to use.
+        mesh_file_name: Name of mesh XML in configuration directory.
+        processes: Number of processes to run.
+        nodes: Number of nodes to run on. Only applicable when running on a multi-node
+            system.
+        cpus_per_process: Number of processing units to use per process. Only
+            applicable when running on a multi-node system.
+        wall_time: Maximum time to allow job to run for. Only applicable when submitting
+            to a job scheduler.
+        neso_outfile: Name of the output file the NESO solver produces, to copy across
+            to the local directory.
+
+    """
+    # If we're running remotely, tell the submission to wait until done
+    if "archer2" in fab.env.remote:
+        fab.update_environment(
+            {
+                "job_dispatch": "cd /work/$project/$project/$username ; sbatch --wait",
+            }
+        )
+
+    # Make a temporary directory to store the configs in - they can be deleted when the
+    # run is done
+    with TemporaryDirectory(dir=config) as temporary_config_dir:
+        config_path = Path(temporary_config_dir)
+
+        # Copy the mesh and conditions files here for FabSIM to read it properly
+        shutil.copyfile(mesh_file, config_path / Path(mesh_file).name)
+        shutil.copyfile("conditions.xml", config_path / "conditions.xml")
+
+        # Run an instance of NESO
+        neso(
+            temporary_config_dir,
+            solver=solver,
+            mesh_file_name=Path(mesh_file).name,
+            processes=processes,
+            nodes=nodes,
+            cpus_per_process=cpus_per_process,
+            wall_time=wall_time,
+        )
+
+    # Retrieve the results
+    fab.fetch_results()
+
+    # Move the results to the local VVUQ directory for further processing
+    local_results_dir = Path(fab.env.job_results_local) / template(
+        fab.env.job_name_template
+    )
+    shutil.copyfile(local_results_dir / neso_outfile, neso_outfile)
+
+
+@fab.task
+@fab.load_plugin_env_vars("FabNESO")
+def neso_vvuq_campaign(  # noqa: PLR0913
+    vvuq_config: str,
+    conditions_template: str = "conditions.template",
+    vvuq_parameters: str = "vvuq_parameters.txt",
+    vvuq_script: str = "FabNESO/run_vvuq_instance.sh",
+    neso_mesh_file: str = "mesh.xml",
+    output_dir: str = "neso_pce/",
+    neso_outfile: str = "Electrostatic2D3V_line_field_evaluations.h5part",
+    solver: str = "Electrostatic2D3V",
+    processes: str | int = 4,
+    nodes: str | int = 1,
+    cpus_per_process: str | int = 1,
+    wall_time: str = "00:15:00",
+    *,
+    run_vvuq_sequential: bool = True,
+) -> None:
+    """
+    Run a VVUQ campaign on a NESO solver.
+
+    Args:
+        vvuq_config: Directory with VVUQ parameters, conditions template file, and mesh
+            configuration.
+
+    Keyword Args:
+        conditions_template: A template NESO conditions file. Parameters that are to be
+            varied are indicated with a $ sign.
+        vvuq_parameters: A file in the vvuq_config directory containing a dict with
+            "parameters" that correspond to the parameters to replace in template file,
+            "output_columns" that are the NESO outputs that VVUQ should read for its
+            analysis
+        neso_mesh_file: Name of mesh XML in vvuq_config directory.
+        output_dir: Directory in which to carry out the VVUQ campaign.
+        neso_outfile: Name of the output file the NESO solver produces, to copy across
+            to the local directory.
+        solver: Which NESO solver to use.
+        processes: Number of processes to run.
+        nodes: Number of nodes to run on. Only applicable when running on a multi-node
+            system.
+        cpus_per_process: Number of processing units to use per process. Only
+            applicable when running on a multi-node system.
+        wall_time: Maximum time to allow job to run for. Only applicable when submitting
+            to a job scheduler.
+        run_vvuq_sequential: Set to False to parallelise. Default is True because this
+            can cause issues if unchecked.
+
+    """
+    path_to_vvuq_configs = Path(fab.find_config_file_path(vvuq_config))
+
+    # Extract the parameters for the conditions template and variations from the config
+    params, vary, outputs_to_read = _parse_vvuq_parameters(
+        path_to_vvuq_configs / vvuq_parameters
+    )
+
+    # Make an encoder
+    encoder = uq.encoders.GenericEncoder(
+        path_to_vvuq_configs / conditions_template,
+        delimiter="$",
+        target_filename="conditions.xml",
+    )
+
+    # Make the VVUQ decoder - this should be more generic!
+    decoder = uq.decoders.HDF5(
+        target_filename=neso_outfile,
+        output_columns=outputs_to_read,
+    )
+
+    # Get the path for this plugin,
+    vvuq_run_path = Path(fab.get_plugin_path("FabNESO")) / vvuq_script
+
+    # The script that runs the individual NESO instances
+    execute_local = uq.actions.ExecuteLocal(
+        f"{vvuq_run_path} {fab.env.remote} {vvuq_config} "
+        f"{path_to_vvuq_configs / neso_mesh_file} {solver} {processes} "
+        f"{nodes} {cpus_per_process} {wall_time} {neso_outfile}"
+    )
+
+    # Put together the actions required for the VVUQ campaign
+    actions = uq.actions.Actions(
+        uq.actions.CreateRunDirectory("/tmp", flatten=True),  # noqa: S108
+        uq.actions.Encode(encoder),
+        execute_local,
+        uq.actions.Decode(decoder),
+    )
+
+    # Make the VVUQ campaign
+    campaign = uq.Campaign(name=output_dir, params=params, actions=actions)
+
+    # Make and attach a sampler to the campaign
+    my_sampler = uq.sampling.PCESampler(vary)
+    campaign.set_sampler(my_sampler)
+
+    # Run VVUQ
+    campaign.execute(sequential=run_vvuq_sequential).collate()
diff --git a/config_files/vvuq_example/conditions.template b/config_files/vvuq_example/conditions.template
new file mode 100644
index 0000000..9f68042
--- /dev/null
+++ b/config_files/vvuq_example/conditions.template
@@ -0,0 +1,77 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+    xsi:noNamespaceSchemaLocation="http://www.nektar.info/schema/nektar.xsd">
+
+    <EXPANSIONS>
+        <E COMPOSITE="C[1]" NUMMODES="2" TYPE="MODIFIED" FIELDS="u" />
+        <E COMPOSITE="C[1]" NUMMODES="2" TYPE="MODIFIED" FIELDS="rho" />
+    </EXPANSIONS>
+
+    <CONDITIONS>
+
+        <SOLVERINFO>
+            <I PROPERTY="EQTYPE" VALUE="PoissonPIC" />
+            <I PROPERTY="Projection" VALUE="Continuous" />
+        </SOLVERINFO>
+
+        <GLOBALSYSSOLNINFO>
+            <V VAR="u">
+            <I PROPERTY="GlobalSysSoln" VALUE="IterativeStaticCond" />
+            <I PROPERTY="IterativeSolverTolerance" VALUE="1e-8"/>
+            </V>
+        </GLOBALSYSSOLNINFO>
+
+        <PARAMETERS>
+            <P> Lambda = 0.0 </P>
+            <P> epsilon = -1.0 </P>
+            <P> num_particles_total = $NUMBER_PARTICLES </P>
+            <P> num_particles_per_cell = -1 </P>
+            <P> particle_time_step = 0.001 </P>
+            <P> particle_num_time_steps = $TIME_STEPS </P>
+            <P> particle_num_write_particle_steps = 0 </P>
+            <P> particle_num_write_field_energy_steps = 20 </P>
+            <P> particle_num_write_field_steps = $WRITE_FIELD_STEPS </P>
+            <P> particle_num_print_steps = 40 </P>
+            <P> particle_distribution_position = 2 </P>
+            <P> particle_initial_velocity = $INIT_VELOCITY </P>
+            <P> particle_charge_density = $CHARGE_DENSITY </P>
+            <P> particle_number_density = $NUMBER_DENSITY </P>
+            <P> line_field_deriv_evaluations_step = $LINE_STEP </P>
+            <P> line_field_deriv_evaluations_numx = $LINE_NUMX </P>
+            <P> line_field_deriv_evaluations_numy = $LINE_NUMY </P>
+        </PARAMETERS>
+
+        <VARIABLES>
+            <V ID="0"> u </V>
+            <V ID="1"> rho </V>
+        </VARIABLES>
+
+        <BOUNDARYREGIONS>
+            <B ID="1"> C[100] </B>
+            <B ID="2"> C[200]  </B>
+            <B ID="3"> C[300] </B>
+            <B ID="4"> C[400] </B>
+        </BOUNDARYREGIONS>
+
+        <BOUNDARYCONDITIONS>
+            <REGION REF="1">
+                <P VAR="u" VALUE="[3]" />
+                <P VAR="rho" VALUE="[3]" />
+            </REGION>
+            <REGION REF="2">
+                <P VAR="u" VALUE="[4]" />
+                <P VAR="rho" VALUE="[4]" />
+            </REGION>
+            <REGION REF="3">
+                <P VAR="u" VALUE="[1]" />
+                <P VAR="rho" VALUE="[1]" />
+            </REGION>
+            <REGION REF="4">
+                <P VAR="u" VALUE="[2]" />
+                <P VAR="rho" VALUE="[2]" />
+            </REGION>
+        </BOUNDARYCONDITIONS>
+
+    </CONDITIONS>
+
+</NEKTAR>
diff --git a/config_files/vvuq_example/mesh.xml b/config_files/vvuq_example/mesh.xml
new file mode 100644
index 0000000..dba3c9f
--- /dev/null
+++ b/config_files/vvuq_example/mesh.xml
@@ -0,0 +1,155 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <GEOMETRY DIM="2" SPACE="2">
+        <VERTEX>
+            <V ID="0">0.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="1">5.00000000e-02 0.00000000e+00 0.00000000e+00</V>
+            <V ID="2">5.00000000e-02 1.00000000e-02 0.00000000e+00</V>
+            <V ID="3">0.00000000e+00 1.00000000e-02 0.00000000e+00</V>
+            <V ID="4">1.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="5">1.00000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="6">1.50000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="7">1.50000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="8">2.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="9">2.00000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="10">2.50000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="11">2.50000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="12">3.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="13">3.00000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="14">3.50000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="15">3.50000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="16">4.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="17">4.00000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="18">4.50000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="19">4.50000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="20">5.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="21">5.00000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="22">5.50000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="23">5.50000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="24">6.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="25">6.00000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="26">6.50000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="27">6.50000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="28">7.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="29">7.00000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="30">7.50000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="31">7.50000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="32">8.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="33">8.00000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="34">8.50000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="35">8.50000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="36">9.00000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="37">9.00000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="38">9.50000000e-01 0.00000000e+00 0.00000000e+00</V>
+            <V ID="39">9.50000000e-01 1.00000000e-02 0.00000000e+00</V>
+            <V ID="40">1.00000000e+00 0.00000000e+00 0.00000000e+00</V>
+            <V ID="41">1.00000000e+00 1.00000000e-02 0.00000000e+00</V>
+        </VERTEX>
+        <EDGE>
+            <E ID="0">0 1</E>
+            <E ID="1">1 2</E>
+            <E ID="2">2 3</E>
+            <E ID="3">3 0</E>
+            <E ID="4">1 4</E>
+            <E ID="5">4 5</E>
+            <E ID="6">5 2</E>
+            <E ID="7">4 6</E>
+            <E ID="8">6 7</E>
+            <E ID="9">7 5</E>
+            <E ID="10">6 8</E>
+            <E ID="11">8 9</E>
+            <E ID="12">9 7</E>
+            <E ID="13">8 10</E>
+            <E ID="14">10 11</E>
+            <E ID="15">11 9</E>
+            <E ID="16">10 12</E>
+            <E ID="17">12 13</E>
+            <E ID="18">13 11</E>
+            <E ID="19">12 14</E>
+            <E ID="20">14 15</E>
+            <E ID="21">15 13</E>
+            <E ID="22">14 16</E>
+            <E ID="23">16 17</E>
+            <E ID="24">17 15</E>
+            <E ID="25">16 18</E>
+            <E ID="26">18 19</E>
+            <E ID="27">19 17</E>
+            <E ID="28">18 20</E>
+            <E ID="29">20 21</E>
+            <E ID="30">21 19</E>
+            <E ID="31">20 22</E>
+            <E ID="32">22 23</E>
+            <E ID="33">23 21</E>
+            <E ID="34">22 24</E>
+            <E ID="35">24 25</E>
+            <E ID="36">25 23</E>
+            <E ID="37">24 26</E>
+            <E ID="38">26 27</E>
+            <E ID="39">27 25</E>
+            <E ID="40">26 28</E>
+            <E ID="41">28 29</E>
+            <E ID="42">29 27</E>
+            <E ID="43">28 30</E>
+            <E ID="44">30 31</E>
+            <E ID="45">31 29</E>
+            <E ID="46">30 32</E>
+            <E ID="47">32 33</E>
+            <E ID="48">33 31</E>
+            <E ID="49">32 34</E>
+            <E ID="50">34 35</E>
+            <E ID="51">35 33</E>
+            <E ID="52">34 36</E>
+            <E ID="53">36 37</E>
+            <E ID="54">37 35</E>
+            <E ID="55">36 38</E>
+            <E ID="56">38 39</E>
+            <E ID="57">39 37</E>
+            <E ID="58">38 40</E>
+            <E ID="59">40 41</E>
+            <E ID="60">41 39</E>
+        </EDGE>
+        <ELEMENT>
+            <Q ID="0">0 1 2 3</Q>
+            <Q ID="1">4 5 6 1</Q>
+            <Q ID="2">7 8 9 5</Q>
+            <Q ID="3">10 11 12 8</Q>
+            <Q ID="4">13 14 15 11</Q>
+            <Q ID="5">16 17 18 14</Q>
+            <Q ID="6">19 20 21 17</Q>
+            <Q ID="7">22 23 24 20</Q>
+            <Q ID="8">25 26 27 23</Q>
+            <Q ID="9">28 29 30 26</Q>
+            <Q ID="10">31 32 33 29</Q>
+            <Q ID="11">34 35 36 32</Q>
+            <Q ID="12">37 38 39 35</Q>
+            <Q ID="13">40 41 42 38</Q>
+            <Q ID="14">43 44 45 41</Q>
+            <Q ID="15">46 47 48 44</Q>
+            <Q ID="16">49 50 51 47</Q>
+            <Q ID="17">52 53 54 50</Q>
+            <Q ID="18">55 56 57 53</Q>
+            <Q ID="19">58 59 60 56</Q>
+        </ELEMENT>
+        <CURVED />
+        <COMPOSITE>
+            <C ID="1"> Q[0-19] </C>
+            <C ID="100"> E[0,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58] </C>
+            <C ID="200"> E[59] </C>
+            <C ID="300"> E[2,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60] </C>
+            <C ID="400"> E[3] </C>
+        </COMPOSITE>
+        <DOMAIN>
+            <D ID="0"> C[1] </D>
+        </DOMAIN>
+    </GEOMETRY>
+    <Metadata>
+        <Provenance>
+            <GitBranch></GitBranch>
+            <GitSHA1>f1598d5e39f175acf388b90df392f76ff29d7f9d</GitSHA1>
+            <Hostname>L0211-XU</Hostname>
+            <NektarVersion>5.2.0</NektarVersion>
+            <Timestamp>23-Nov-2022 10:54:06</Timestamp>
+        </Provenance>
+        <NekMeshCommandLine>-v -m peralign:surf1=100:surf2=300:dir=y -m peralign:surf1=200:surf2=400:dir=x periodic_structured_cartesian_rectangle.msh periodic_structured_cartesian_rectangle.xml:xml:uncompress </NekMeshCommandLine>
+    </Metadata>
+</NEKTAR>
diff --git a/config_files/vvuq_example/vvuq_parameters.txt b/config_files/vvuq_example/vvuq_parameters.txt
new file mode 100644
index 0000000..8ca9311
--- /dev/null
+++ b/config_files/vvuq_example/vvuq_parameters.txt
@@ -0,0 +1,17 @@
+{
+ "parameters":
+ {"NUMBER_PARTICLES" : [40000],
+  "INIT_VELOCITY" : [1.0,0.9,1.1],
+  "CHARGE_DENSITY" : [105,101,109],
+  "NUMBER_DENSITY" : [105],
+  "TIME_STEPS" : [180],
+  "WRITE_FIELD_STEPS": [100],
+  "LINE_STEP": [20],
+  "LINE_NUMX": [100],
+  "LINE_NUMY": [1],
+ },
+ "output_columns": [
+  "Step#160/FIELD_EVALUATION_0",
+  "Step#160/x",
+ ]
+}
diff --git a/requirements.txt b/requirements.txt
index 167df30..eae1bf5 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,2 +1,4 @@
 h5py>=3.10.0
 PyVBMC>=1.0.1
+easyvvuq>=1.2
+chaospy>=4.3.2