diff --git a/docs/examples/simulator/plot_run_ravaflow24.py b/docs/examples/simulator/plot_run_ravaflow24.py
new file mode 100644
index 0000000000000000000000000000000000000000..96a2cbc0e87964bdd63329ccc5d1afeb007c998a
--- /dev/null
+++ b/docs/examples/simulator/plot_run_ravaflow24.py
@@ -0,0 +1,159 @@
+"""
+
+RunSimulator: Ravaflow24 Mixture Model
+======================================
+"""
+
+# %% md
+# 
+# This example shows how to run multiple Ravaflow24Mixture simulations in serial
+# and parallelly using :class:`.RunSimulator`.
+#
+
+# %% md
+# 
+# First, we import the class :class:`.Ravaflow24Mixture`.
+
+from psimpy.simulator import Ravaflow24Mixture
+
+# %% md
+# 
+# To create an instance of this class, we must specify the parameter ``dir_sim``.
+# It represents the directory in which output files generated by r.avaflow will be
+# saved.
+
+import os
+
+# Here we create a folder called `temp1_run_Ravaflow24Mixture_example` to save
+# output files generated by r.avaflow
+cwd = os.getcwd()
+if not os.path.exists('temp1_run_Ravaflow24Mixture_example'):
+    os.mkdir('temp1_run_Ravaflow24Mixture_example')
+dir_sim = os.path.join(cwd, 'temp1_run_Ravaflow24Mixture_example')
+
+# %% md
+# 
+# Now, we can create an instance of :class:`.Ravaflow24Mixture` by given
+# ``dir_sim``. To reduce simulation time, we set ``time_end`` to :math:`50`.
+# Other parameters are set to their default values.
+
+voellmy_model = Ravaflow24Mixture(dir_sim, time_end=50)
+
+# %% md
+# 
+# In this example, we define a `simulator` that takes required inputs and returns
+# the overall impact area based on ``voellmy_model``.
+
+import numpy as np
+
+def simulator(prefix, elevation, hrelease, basal_friction, turbulent_friction, EPSG):
+    """Preprocess required inputs, run simulation, and return output as a numpy array."""
+    grass_location, sh_file = voellmy_model.preprocess(prefix=prefix, elevation=elevation,
+        hrelease=hrelease, basal_friction=basal_friction, turbulent_friction=turbulent_friction, EPSG=EPSG)
+    voellmy_model.run(grass_location, sh_file) # running this line, r.avaflow will write outputs to dir_sim
+    impact_area = voellmy_model.extract_impact_area(prefix)
+    
+    return np.array([impact_area]) # define dir_out and set save_out to True, returned numpy array will be saved to dir_out
+
+
+# %% md
+# 
+# We are going to run multiple simulations at different values of ``basal_friction``
+# and ``turbulent_friction``.
+
+import itertools
+
+var_inp_parameter = ['basal_friction', 'turbulent_friction']
+basal_friction = [20, 30] 
+turbulent_friction = [3, 4]
+var_samples = np.array(
+    [x for x in itertools.product(basal_friction, turbulent_friction)])
+
+
+# %% md
+#
+# Other parameters are fixed. 
+
+dir_data = os.path.abspath('../../../tests/data/')
+elevation = os.path.join(dir_data, 'synthetic_topo.tif')
+hrelease = os.path.join(dir_data, 'synthetic_rel.tif')
+fix_inp = {'elevation': elevation, 'hrelease': hrelease, 'EPSG': '2326'}
+
+# %% md
+# 
+# The parameter ``prefix`` of the `simulator` is special. It is not involved in
+# the computational model, but only used to name output files generated by r.avaflow.
+# Usually, we need a unique value of ``prefix`` for each simulation to avoid conflict.
+# Such parameter is not defined in ``var_inp_parameter`` or ``fix_inp`` of
+# :class:`.RunSimulator`. Instead, we use a seperate parameter, called ``o_parameter``
+# for this purpose.
+
+o_parameter = 'prefix'
+
+
+# %% md
+# 
+# We may want to save outputs returned by `simulator` at each simulation for
+# later inspection or processing. In that case, we need to define ``dir_out`` and
+# set ``save_out`` as `True`.
+
+import os
+
+# Here we create a folder called `temp2_run_Ravaflow24Mixture_example` to save outputs
+# returned at each simulation.
+cwd = os.getcwd()
+if not os.path.exists('temp2_run_Ravaflow24Mixture_example'):
+    os.mkdir('temp2_run_Ravaflow24Mixture_example')
+dir_out = os.path.join(cwd, 'temp2_run_Ravaflow24Mixture_example')
+
+# %% md
+#
+# .. warning:: ``dir_sim`` and ``dir_out`` are two different directories for
+#    different purposes. ``dir_sim`` stores outputs generated by the function body
+#    of `simulator`, usually a third-party software or solver (here r.avaflow).
+#    ``dir_out`` stores returned numpy array of `simulator`, which include outputs
+#    of interest (here impact area). They can be the same directory if file names 
+#    have no conflict. We recommend to keep them seperate. In addition, if the function
+#    body of `simulator` does not save files to disk, ``dir_sim`` is not needed. Our
+#    :class:`.MassPointModel` is an example. Similarly, if we do not want to save
+#    returned numpy array of `simulator`, ``dir_out`` is not needed.
+
+# %% md
+#
+# Now we can define an object of :class:`.RunSimulator` by
+
+from psimpy.simulator import RunSimulator
+
+run_simulator = RunSimulator(simulator, var_inp_parameter, fix_inp, o_parameter,
+    dir_out, save_out=True)
+
+
+# %% md
+# 
+# Define prefixes which will be used to name files generated by each r.avaflow
+# simulation and each returned numpy array of `simulator`.
+
+serial_prefixes = ["serial"+str(i) for i in range(len(var_samples))]
+parallel_prefixes = ["parallel"+str(i) for i in range(len(var_samples))]
+
+import time
+
+start = time.time()
+run_simulator.serial_run(var_samples=var_samples, prefixes=serial_prefixes)
+serial_time = time.time() - start
+serial_output = run_simulator.outputs
+print(f"serial_output: {serial_output}")
+
+start = time.time()
+run_simulator.parallel_run(var_samples, prefixes=parallel_prefixes,
+    max_workers=2, append=True)
+parallel_time = time.time() - start
+    
+parallel_output = run_simulator.outputs[len(var_samples):]
+print(f"parallel_output: {parallel_output}")
+
+    
+print("Serial run time: ", serial_time)
+print("Parallel run time: ", parallel_time)
+
+