diff --git a/docs/source/api.rst b/docs/source/api.rst
index 04faab1..1916abd 100644
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -10,6 +10,18 @@ EchemSolver classes
 .. autoclass:: echemfem.TransientEchemSolver
    :members:
 
+FlowSolver classes
+-------------------
+
+.. autoclass:: echemfem.FlowSolver
+   :members:
+
+.. autoclass:: echemfem.NavierStokesFlowSolver
+   :members:
+
+.. autoclass:: echemfem.NavierStokesBrinkmanFlowSolver
+   :members:
+
 Utility meshes
 --------------
 These are boundary layer meshes, which are modified versions of :func:`firedrake.utility_meshes.IntervalMesh` and :func:`firedrake.utility_meshes.RectangleMesh`.
diff --git a/docs/source/user_guide/fluid_solver.rst b/docs/source/user_guide/fluid_solver.rst
new file mode 100644
index 0000000..84205d3
--- /dev/null
+++ b/docs/source/user_guide/fluid_solver.rst
@@ -0,0 +1,65 @@
+Fluid Flow Solvers
+===================
+
+For convenience, we provide a simple implementation of fluid flow equations in
+:class:`echemfem.FlowSolver`, which can be used to obtain a velocity field for
+the :class:`echemfem.EchemSolver`.
+
+Navier-Stokes Solver
+--------------------
+
+The :class:`echemfem.NavierStokesFlowSolver` class contains an implementation
+for the incompressible Navier-Stokes equations. Both nondimensional and
+dimensional equations are available. Physical parameters are passed as a
+:py:class:`dict` in the ``fluid_params`` argument.
+
+To use the nondimensional version, the user must pass the following key:
+
+* ``"Reynolds number"``
+
+For the dimensional version, the user must pass the following keys:
+
+* ``"density"``
+
+* ``"dynamic viscosity"`` or ``"kinematic viscosity"``
+
+Navier-Stokes-Brinkman Solver
+-----------------------------
+
+The :class:`echemfem.NavierStokesBrinkmanFlowSolver` class contains an
+implementation for the incompressible Navier-Stokes-Brinkman equations. Both
+nondimensional and dimensional equations are available. In addition to the
+parameters provided for Navier-Stokes, the physical parameters below must be
+provided.
+
+To use the nondimensional version, the user must also pass the following key:
+
+* ``"Darcy number"``
+
+For the dimensional version, the user must also pass the following keys:
+
+* ``"permeability"`` or ``"inverse permeability"``
+
+* Optional: ``"effective kinematic viscosity"``
+
+Boundary conditions
+-------------------
+
+A :py:class:`dict` containing boundary markers is passed when creating a
+``FlowSolver`` object. Below are the different options for ``boundary_markers``
+keys. The velocity is denoted as :math:`\mathbf u` and pressure as :math:`p`.
+
+* ``"no slip"``: Sets velocity to zero, :math:`\mathbf u = 0`.
+
+* ``"inlet velocity"``: Sets inlet velocity, :math:`\mathbf u = \mathbf
+  u_\mathrm{in}`, which is passed in ``fluid_params`` as ``"inlet velocity"``.
+
+* ``"outlet velocity"``: Sets outlet velocity, :math:`\mathbf u = \mathbf
+  u_\mathrm{out}`, which is passed in ``fluid_params`` as ``"outlet
+  velocity"``.
+
+* ``"inlet pressure"``: Sets inlet velocity, :math:`p = p_\mathrm{in}`, which
+  is passed in ``fluid_params`` as ``"inlet pressure"``.
+
+* ``"outlet pressure"``: Sets outlet pressure, :math:`p = p_\mathrm{out}`,
+  which is passed in ``fluid_params`` as ``"outlet pressure"``.
diff --git a/docs/source/user_guide/index.rst b/docs/source/user_guide/index.rst
index ddaa04c..545dc3e 100644
--- a/docs/source/user_guide/index.rst
+++ b/docs/source/user_guide/index.rst
@@ -59,3 +59,11 @@ For transient cases, a temporal grid defined as a  ``list`` or
    import numpy as np
    times = np.linspace(0, 11, 1)  # 10 timesteps of size 0.1
    solver.solve(times)
+
+To generate non-trivial velocity fields, some fluid flow solvers are available:
+
+.. toctree::
+   :maxdepth: 2
+
+   fluid_solver
+
diff --git a/echemfem/__init__.py b/echemfem/__init__.py
index fad0b4c..a9561f5 100644
--- a/echemfem/__init__.py
+++ b/echemfem/__init__.py
@@ -1,4 +1,4 @@
 from echemfem.cylindricalmeasure import CylindricalMeasure
 from echemfem.solver import EchemSolver, TransientEchemSolver
-from echemfem.flow_solver import NavierStokesFlowSolver, NavierStokesBrinkmanFlowSolver
+from echemfem.flow_solver import FlowSolver, NavierStokesFlowSolver, NavierStokesBrinkmanFlowSolver
 from echemfem.utility_meshes import IntervalBoundaryLayerMesh, RectangleBoundaryLayerMesh
diff --git a/echemfem/flow_solver.py b/echemfem/flow_solver.py
index 39cf965..f62e51b 100644
--- a/echemfem/flow_solver.py
+++ b/echemfem/flow_solver.py
@@ -7,18 +7,19 @@ class FlowSolver(ABC):
     """Base class for a flow solver.
 
     This class is used to create solvers for fluid flow decoupled from the chemistry
+
     Attributes:
         mesh (:class:`firedrake.mesh.MeshGeometry`): Mesh object from firedrake
-        flow_params (dict): Dictionary containing physical parameters 
+        fluid_params (dict): Dictionary containing physical parameters 
         boundary_markers (dict): Dictionary where the keys are :py:class:`str:`
-        representing the type of boundary condition, and the values are
-        :py:class:`tuple` containing the boundary indices. For example :
+            representing the type of boundary condition, and the values are
+            :py:class:`tuple` containing the boundary indices. For example :
 
-        .. code-block::
+            .. code-block::
 
-            self.boundary_markers = {"inlet velocity": (1,)}
+                boundary_markers = {"inlet velocity": (1,)}
 
-        sets the boundary condition ``"inlet velocity"`` on boundary 1.
+            sets the boundary condition ``"inlet velocity"`` on boundary 1.
     """
 
     def __init__(self, mesh, fluid_params, boundary_markers):
@@ -174,6 +175,7 @@ def setup_problem(self):
 
     def setup_solver(self, ksp_solver="lu"):
         """Optional PCD preconditioner for nondimensional Navier-Stokes
+
         Args:
             ksp_solver (str): ``"lu"`` or ``"pcd"``
         """
@@ -221,8 +223,9 @@ class NavierStokesBrinkmanFlowSolver(FlowSolver):
 
     """Incompressible Navier-Stokes-Brinkman solver
 
-       For dimensionless form, pass Reynolds number to fluid_params.
-       For dimensional form, pass density and kinematic viscosity.
+       For dimensionless form, pass Reynolds and Darcy numbers to fluid_params.
+       For dimensional form, pass density, permeability, and kinematic or
+       dynamic viscosity.
     """
     def setup_problem(self):
         u = self.u