diff --git a/docs/source_docs/refs.bib b/docs/source_docs/refs.bib
index 63e613739ae4158edfe63f4bc588a6aebcf36ed3..6fa2b3d9b9ff4e568e7353eefa35307e30818e65 100644
--- a/docs/source_docs/refs.bib
+++ b/docs/source_docs/refs.bib
@@ -291,7 +291,6 @@ url = {https://www.sciencedirect.com/science/article/pii/S0304885312002132},
 }
 
 
-
 @article{garg12,
   title={Open-source {MFIX-DEM} software for gas-solids flows: {Part I-Verification studies}},
   author={Garg, Rahul and Galvin, Janine and Li, Tingwen and Pannala, Sreekanth},
@@ -302,6 +301,18 @@ url = {https://www.sciencedirect.com/science/article/pii/S0304885312002132},
   publisher={Elsevier}
 }
 
+
+@article{gartling90,
+  author = {Gartling, David K.},
+  title = {A test problem for outflow boundary conditions—flow over a backward-facing step},
+  journal = {International Journal for Numerical Methods in Fluids},
+  volume = {11},
+  number = {7},
+  pages = {953-967},
+  doi = {https://doi.org/10.1002/fld.1650110704},
+  year = {1990}
+}
+
 @article{gibilaro07,
   author = {Gibilaro L.G., Gallucci K., Felice R. Di. and Pagliai P.},
   title = {On the apparent viscosity of a fluidized bed},
@@ -363,6 +374,14 @@ doi = {https://doi.org/10.1016/j.jcp.2022.111305},
   pages = {843-855}
 }
 
+@article{johnson99,
+  title={Flow past a sphere up to a Reynolds number of 300},
+  volume={378}, DOI={10.1017/S0022112098003206},
+  journal={Journal of Fluid Mechanics},
+  author={JOHNSON, T. A. and PATEL, V. C.},
+  year={1999}, pages={19–70}
+}
+
 @article{krieger59,
   author = {Krieger, I.M. and Dougherty, T.J. },
   title = {A mechanism for non-newtonian flow in suspensions of rigid spheres},
@@ -595,6 +614,17 @@ doi = {https://doi.org/10.1016/j.jcp.2022.111305},
   year = {1988}
 }
 
+@article{taneda56,
+  author = {Taneda ,Sadatoshi},
+  title = {Experimental Investigation of the Wake behind a Sphere at Low Reynolds Numbers},
+  journal = {Journal of the Physical Society of Japan},
+  volume = {11},
+  number = {10},
+  pages = {1104-1108},
+  year = {1956},
+  doi = {10.1143/JPSJ.11.1104}
+}
+
 @article{tang15,
    author = {Tang, Y. and Peters, E. A. J. F. and Kuipers, J. A. M. and Kriebitzsch, S. H. L. and van der Hoef, M. A.},
    title = {A New Drag Correlation from Fully Resolved Simulations of Flow Past Monodisperse Static Arrays of Spheres},
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-richardson_L1.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-richardson_L1.png
new file mode 100644
index 0000000000000000000000000000000000000000..0883fbfff179e19f03577709225e42bb0e11b04d
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-richardson_L1.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-richardson_Lu.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-richardson_Lu.png
new file mode 100644
index 0000000000000000000000000000000000000000..30bb5a5dc7f9c728c6ed076f27f38d1cb940ff84
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-richardson_Lu.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-streamfunction_contours-no-fill.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-streamfunction_contours-no-fill.png
new file mode 100644
index 0000000000000000000000000000000000000000..a314569a61eaf53c3f28f030181aa9df93af5c8f
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-streamfunction_contours-no-fill.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-streamfunction_contours.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-streamfunction_contours.png
new file mode 100644
index 0000000000000000000000000000000000000000..90444efb33d213c8de9f4c45e55ba3ff442c0ccd
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-streamfunction_contours.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-u_velocity_all_cases.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-u_velocity_all_cases.png
new file mode 100644
index 0000000000000000000000000000000000000000..8df979e57f4f46250c79e1ebc5b6b22b36929e3e
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-u_velocity_all_cases.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-v_velocity_all_cases.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-v_velocity_all_cases.png
new file mode 100644
index 0000000000000000000000000000000000000000..1aa814741faa0c1d5d32c5a268f94ea4f2230e05
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-v_velocity_all_cases.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-vorticity_contours-no-fill.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-vorticity_contours-no-fill.png
new file mode 100644
index 0000000000000000000000000000000000000000..93d5e47d2c7ea6d88291d860db9733e5d4ffbacb
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-vorticity_contours-no-fill.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-vorticity_contours.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-vorticity_contours.png
new file mode 100644
index 0000000000000000000000000000000000000000..9de77614a9f1d2c3a465ce53664d56aadd4408d6
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-vorticity_contours.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-wall_pressure_profiles.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-wall_pressure_profiles.png
new file mode 100644
index 0000000000000000000000000000000000000000..720e0a21de22dab809a740782eb3b84ae025c935
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-backward-facing-step-wall_pressure_profiles.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/images/steady-flow-past-sphere-midplane_contour.png b/docs/source_docs/tutorials/incompressible-fluid/images/steady-flow-past-sphere-midplane_contour.png
new file mode 100644
index 0000000000000000000000000000000000000000..c4a302e0d1813f301deaa884c3ac5e4f58fd07b4
Binary files /dev/null and b/docs/source_docs/tutorials/incompressible-fluid/images/steady-flow-past-sphere-midplane_contour.png differ
diff --git a/docs/source_docs/tutorials/incompressible-fluid/steady-backward-facing-step.rst b/docs/source_docs/tutorials/incompressible-fluid/steady-backward-facing-step.rst
new file mode 100644
index 0000000000000000000000000000000000000000..503571f9fac1531bd751d0ba27a084d264335f06
--- /dev/null
+++ b/docs/source_docs/tutorials/incompressible-fluid/steady-backward-facing-step.rst
@@ -0,0 +1,421 @@
+.. _TutorialSteadyBackwardFacingStep:
+
+Steady flow over a backward-facing step
+=======================================
+
+This tutorial sets up incompressible flow over a backward-facing step based
+on the benchmark problem of Gartling :cite:p:`gartling90`. This is a steady,
+Reynolds-800 channel-flow example for examining recirculation downstream of
+a sudden expansion and the effect of the pressure-outlet boundary condition.
+
+Features
+--------
+
+* Incompressible fluid with no particles.
+* Steady-state advance.
+* Backward-facing step represented by the inlet velocity profile.
+* Parabolic velocity inlet, pressure outlet, no-slip walls, and a periodic
+  spanwise direction.
+* Plotfile output for velocity and pressure.
+* Comparison plots for velocity profiles, wall pressure, streamfunction, and
+  vorticity, plus a grid-convergence study using Richardson extrapolation.
+
+Case description
+----------------
+
+The computational domain begins at the step and extends downstream from
+:math:`x=0` to :math:`x=30`. The expanded channel occupies
+:math:`-0.5 \le y \le 0.5`. Flow enters through the upper half of the
+low-:math:`x` boundary, :math:`0 < y \le 0.5`, and the lower half of that
+boundary has zero inlet velocity to represent the step face. The low- and
+high-:math:`y` boundaries are no-slip walls. To make the case quasi-two
+dimensional, the domain is periodic in :math:`z`.
+
+To bound the cost of the spanwise direction as the streamwise and
+cross-channel resolution is refined, the :math:`z`-domain is sized so that
+the mesh cells remain cubic, :math:`\Delta x = \Delta y = \Delta z`. At the
+base resolution the number of cells across :math:`z` equals the number of
+cells across :math:`y`, and the :math:`z`-domain spans the full channel
+height, :math:`-0.5 \le z \le 0.5`. For the finer meshes used in the
+grid-convergence study below, the number of cells across :math:`z` is capped
+at 32 and the :math:`z`-domain length is reduced accordingly so that
+:math:`\Delta z` continues to match :math:`\Delta x` and :math:`\Delta y`.
+Applying this rule gives the following domain sizes:
+
+.. list-table:: Spanwise (z) domain sizing across mesh resolutions
+   :header-rows: 1
+
+   * - Mesh
+     - :math:`n_y`
+     - :math:`n_z`
+     - :math:`\Delta x = \Delta y = \Delta z`
+     - :math:`z`-domain length
+   * - 480x16
+     - 16
+     - 16
+     - 0.0625
+     - 1.000
+   * - 960x32
+     - 32
+     - 32
+     - 0.03125
+     - 1.000
+   * - 1920x64
+     - 64
+     - 32
+     - 0.015625
+     - 0.500
+   * - 3840x128
+     - 128
+     - 32
+     - 0.0078125
+     - 0.250
+   * - 7680x256
+     - 256
+     - 32
+     - 0.00390625
+     - 0.125
+
+The tutorial uses
+
+.. math::
+
+   \rho = 1,\qquad \mu = 0.00125,\qquad H = 1,
+
+where :math:`H` is the expanded channel height. The inlet velocity is
+
+.. math::
+
+   u(x=0,y) =
+   \begin{cases}
+   24y(0.5-y), & y > 0,\\
+   0,          & y \le 0.
+   \end{cases}
+
+The mean velocity over the inlet opening is :math:`\overline{u}=1`, giving
+
+.. math::
+
+   Re = \frac{\rho \overline{u} H}{\mu} = 800.
+
+At this Reynolds number the solution contains a primary recirculation region
+behind the lower step and a weaker upper-wall recirculation region farther
+downstream.
+
+
+Important input sections
+------------------------
+
+The mesh is a long channel with a periodic :math:`z` direction:
+
+.. code-block:: text
+
+   geometry.is_periodic =  0    0     1
+   geometry.prob_lo     =  0.  -0.5  -0.5
+   geometry.prob_hi     = 30.   0.5   0.5
+
+   amr.n_cell = 480  16  16
+
+No embedded-boundary geometry is needed for this setup. The step is imposed
+through the inlet condition:
+
+.. code-block:: text
+
+   mfix.geometry = None
+
+The low- and high-:math:`y` faces are no-slip walls. The low-:math:`x` face is
+a mass inflow with the Gartling parabolic profile :cite:p:`gartling90` over the
+upper half of the opening, and the high-:math:`x` face is a pressure outlet:
+
+.. code-block:: text
+
+   bc.regions = bottom_wall top_wall inflow outflow
+
+   bc.bottom_wall = no-slip
+   bc.top_wall = no-slip
+
+   bc.inflow = mi
+   bc.inflow.fluid.density = 1.0
+   bc.inflow.fluid.inflow_type = velocity
+   bc.inflow.fluid.velocity = "if(y>0., 24.*y*(0.5 - y), 0.)"
+
+   bc.outflow = po
+   bc.outflow.pressure = 0.0
+
+The fluid viscosity sets the Reynolds number:
+
+.. code-block:: text
+
+   fluid.viscosity.molecular.model = constant
+   fluid.viscosity.molecular.constant = 0.00125
+
+The case is run with the pseudo steady-state time advance with a tolerance
+of :math:`1.0e-6` and 250,000 maximum iterations.
+
+.. code-block:: text
+
+   mfix.steady_state = 1
+
+   mfix.steady_state_tol = 1.0e-6
+   mfix.steady_state_maxiter = 250000
+
+   mfix.dt_max = 1.e-3
+
+For the grid convergence study, the tolerance is reduced to :math:`1e-8`,
+and the maximum iterations is increased to 500,000.
+
+.. note::
+
+   The complete input file is located in the MFIX-Exa source directory as
+   ``tutorials/incompressible-fluid/inputs.steady-backward-facing-step``
+
+
+Post-Processing
+---------------
+
+The run writes an AMReX plotfile named ``plt*`` every 5000 iterations into the
+``io`` subdirectory. Each plot file saves the fluid:
+
+* velocity field, ``vel_g``, and
+* pressure, ``p_g``.
+
+The plotfile data can be viewed in ParaView, VisIt, or any other
+AMReX-compatible post-processing tool. The streamfunction and vorticity
+contours below were generated in a Jupyter notebook using the yt reader,
+and show the lower-wall separation bubble and the weaker upper-wall
+recirculation region.
+
+.. _fig_steady_backward_facing_step_streamfunction:
+
+.. figure:: ./images/steady-backward-facing-step-streamfunction_contours.png
+   :width: 100%
+   :align: center
+   :alt: Streamfunction contours for steady flow over a backward-facing step.
+
+   Streamfunction contours at steady state.
+
+
+.. _fig_steady_backward_facing_step_vorticity:
+
+.. figure:: ./images/steady-backward-facing-step-vorticity_contours.png
+   :width: 100%
+   :align: center
+   :alt: Vorticity contours for steady flow over a backward-facing step.
+
+   Vorticity contours at steady state.
+
+
+The horizontal velocity profiles at :math:`x=7` and :math:`x=15` can be used
+to compare the MFIX-Exa solution with Gartling's tabulated results :cite:p:`gartling90`.
+The figure below includes several uniformly refined meshes, labeled by the number of
+cells in the :math:`x` and :math:`y` directions.
+
+.. _fig_steady_backward_facing_step_u_profiles:
+
+.. figure:: ./images/steady-backward-facing-step-u_velocity_all_cases.png
+   :width: 70%
+   :align: center
+   :alt: Horizontal velocity profiles at x equals 7 and x equals 15 compared with Gartling 1990 data.
+
+   Horizontal velocity profiles compared with Gartling :cite:p:`gartling90`.
+
+The horizontal velocity is well predicted even on the coarsest mesh: all five
+grids agree closely with Gartling's data :cite:p:`gartling90` at both :math:`x=7` and
+:math:`x=15`, with only modest, monotonically shrinking differences as the
+mesh is refined.
+
+The vertical velocity tells a different story. At :math:`x=15`, downstream of
+reattachment, :math:`v` is small in magnitude and every mesh -- including the
+coarsest -- agrees well with Gartling's data :cite:p:`gartling90`. At :math:`x=7`,
+however, which sits inside the primary recirculation region, the two coarsest
+meshes (480x16 and 960x32) fail to capture the profile. The solution only becomes
+qualitatively correct at 1920x64.
+
+.. _fig_steady_backward_facing_step_v_profiles:
+
+.. figure:: ./images/steady-backward-facing-step-v_velocity_all_cases.png
+   :width: 70%
+   :align: center
+   :alt: Vertical velocity profiles at x equals 7 and x equals 15 compared with Gartling 1990 data.
+
+   Vertical velocity profiles compared with Gartling :cite:p:`gartling90`.
+
+This contrast is expected rather than concerning: :math:`u` is the dominant,
+streamwise velocity component and is comparatively insensitive to resolution,
+while :math:`v` at :math:`x=7` is a small-magnitude, secondary quantity
+generated entirely by the curvature of the shear layer bounding the
+recirculation zone -- the same near-wall region identified in
+`Grid Convergence`_ as needing the finest resolution to correctly locate the
+separation and reattachment points. A grid coarse enough to under-resolve
+that shear layer can still reproduce :math:`u` reasonably well while getting
+:math:`v` qualitatively wrong.
+
+The pressure along the upper and lower walls is another useful comparison
+quantity. The plot below shows the pressure profiles from the finest mesh.
+
+.. _fig_steady_backward_facing_step_wall_pressure:
+
+.. figure:: ./images/steady-backward-facing-step-wall_pressure_profiles.png
+   :width: 70%
+   :align: center
+   :alt: Pressure profiles along the upper and lower channel walls.
+
+   Pressure profiles along the upper and lower channel walls.
+
+
+Grid Convergence
+----------------
+
+To confirm that the reported comparisons are grid converged, the case was run
+on five uniformly refined meshes, labeled by the number of cells in the
+:math:`x` and :math:`y` directions (480x16, 960x32, 1920x64, 3840x128, and
+7680x256). For each mesh, the lower- and upper-wall recirculation zones are
+characterized the same way Gartling :cite:p:`gartling90` does in Tables II and III:
+
+* the vortex center is located as the extremum of the computed streamfunction
+  within each recirculation zone,
+* the streamfunction and vorticity are evaluated at that location, and
+* the separation and reattachment points are located from the zero-crossings
+  of the wall shear stress, :math:`\tau = \mu\, \partial u/\partial y`, using
+  linear interpolation between grid points, following the same method
+  Gartling describes.
+
+Lower-wall eddy
+~~~~~~~~~~~~~~~
+
+.. list-table:: Lower-wall eddy vs. mesh resolution (cf. Gartling :cite:p:`gartling90`, Table II)
+   :header-rows: 1
+
+   * - Mesh
+     - Vortex center :math:`(x, y)`
+     - :math:`\psi`
+     - :math:`\omega`
+     - :math:`L_1`
+   * - 480x16
+     - (2.9062, -0.2188)
+     - -0.0299
+     - -2.2402
+     - 4.9546
+   * - 960x32
+     - (3.2344, -0.2031)
+     - -0.0331
+     - -2.2513
+     - 5.8432
+   * - 1920x64
+     - (3.4453, -0.2109)
+     - -0.0339
+     - -2.2398
+     - 6.0420
+   * - 3840x128
+     - (3.4180, -0.2070)
+     - -0.0341
+     - -2.2488
+     - 6.0869
+   * - 7680x256
+     - (3.4043, -0.2051)
+     - -0.0342
+     - -2.2595
+     - 6.0951
+   * - Gartling (1990)
+     - (3.35, -0.20)
+     - -0.0342
+     - -2.283
+     - 6.10
+
+Richardson extrapolation on the reattachment length :math:`L_1` gives an
+observed order of accuracy of :math:`p \approx 2.1`-:math:`2.5` across the
+three available triplets of meshes, consistent with the nominal second-order
+accuracy of the numerical scheme. Using the finest triplet (1920x64, 3840x128,
+7680x256), :math:`p \approx 2.45` and the extrapolated value is
+:math:`L_1(\Delta x \to 0) \approx 6.097`, matching Gartling's :cite:p:`gartling90`
+reported value of 6.10 to within 0.05%.
+
+.. _fig_steady_backward_facing_step_richardson_L1:
+
+.. figure:: ./images/steady-backward-facing-step-richardson_L1.png
+   :width: 90%
+   :align: center
+   :alt: Richardson extrapolation of the lower-wall reattachment length versus grid spacing.
+
+   Richardson extrapolation of the lower-wall reattachment length :math:`L_1`.
+
+Upper-wall eddy
+~~~~~~~~~~~~~~~
+
+.. list-table:: Upper-wall eddy vs. mesh resolution (cf. Gartling :cite:p:`gartling90`, Table III)
+   :header-rows: 1
+
+   * - Mesh
+     - Vortex center :math:`(x, y)`
+     - :math:`\psi`
+     - :math:`\omega`
+     - Separation :math:`(x, y)`
+     - Reattachment :math:`(x, y)`
+     - :math:`L_u`
+   * - 480x16
+     - (6.9688, 0.2812)
+     - 0.5105
+     - 1.6707
+     - (3.4912, 0.50)
+     - (9.7500, 0.50)
+     - 6.2588
+   * - 960x32
+     - (7.2969, 0.2969)
+     - 0.5072
+     - 1.3218
+     - (4.4186, 0.50)
+     - (10.5215, 0.50)
+     - 6.1029
+   * - 1920x64
+     - (7.4609, 0.3203)
+     - 0.5067
+     - 1.0892
+     - (4.7017, 0.50)
+     - (10.5332, 0.50)
+     - 5.8315
+   * - 3840x128
+     - (7.4492, 0.3164)
+     - 0.5066
+     - 1.1344
+     - (4.7939, 0.50)
+     - (10.5113, 0.50)
+     - 5.7174
+   * - 7680x256
+     - (7.4473, 0.3145)
+     - 0.5065
+     - 1.1584
+     - (4.8268, 0.50)
+     - (10.4965, 0.50)
+     - 5.6697
+   * - Gartling (1990)
+     - (7.40, 0.30)
+     - 0.5064
+     - 1.322
+     - (4.85, 0.50)
+     - (10.48, 0.50)
+     - 5.63
+
+Richardson extrapolation on the upper-wall recirculation length :math:`L_u`
+is more sensitive to mesh resolution than :math:`L_1`: the coarsest triplet
+(480x16, 960x32, 1920x64) gives an unphysical negative order of accuracy,
+which is not unexpected since the upper-wall separation point is a gradual,
+adverse-pressure-gradient-driven separation rather than the fixed, geometric
+separation at the step corner, and is correspondingly harder to resolve on
+coarse meshes. The two finer triplets are consistent with one another,
+:math:`p \approx 1.25`-:math:`1.26`, giving an extrapolated value of
+:math:`L_u(\Delta x \to 0) \approx 5.635`, matching Gartling's :cite:p:`gartling90`
+reported value of 5.63 to within 0.1%.
+
+.. _fig_steady_backward_facing_step_richardson_Lu:
+
+.. figure:: ./images/steady-backward-facing-step-richardson_Lu.png
+   :width: 90%
+   :align: center
+   :alt: Richardson extrapolation of the upper-wall recirculation length versus grid spacing.
+
+   Richardson extrapolation of the upper-wall recirculation length :math:`L_u`.
+
+Taken together, the vortex-center locations and strengths agree with Gartling's
+finest-mesh results :cite:p:`gartling90` to within 1-2% even at moderate resolution,
+and the Richardson-extrapolated separation/reattachment lengths for both
+recirculation zones agree with his reported values to within 0.1%.
diff --git a/docs/source_docs/tutorials/incompressible-fluid/steady-flow-around-cylinder.rst b/docs/source_docs/tutorials/incompressible-fluid/steady-flow-around-cylinder.rst
index 1d331d0baae0d8f3152df3c76f656ad1f7e2e460..389a71131131544cebc516e0c928a79c29f20359 100644
--- a/docs/source_docs/tutorials/incompressible-fluid/steady-flow-around-cylinder.rst
+++ b/docs/source_docs/tutorials/incompressible-fluid/steady-flow-around-cylinder.rst
@@ -1,6 +1,6 @@
 .. _TutorialSteadyFlowAroundCylinder:
 
-Steady Flow Around a Cylinder
+Steady flow around a cylinder
 =============================
 
 This tutorial sets up fluid flow around a fixed circular cylinder based
diff --git a/docs/source_docs/tutorials/incompressible-fluid/steady-flow-around-sphere.rst b/docs/source_docs/tutorials/incompressible-fluid/steady-flow-around-sphere.rst
new file mode 100644
index 0000000000000000000000000000000000000000..1e766db1cf14629c77e64f0c5a704f0ddba509ac
--- /dev/null
+++ b/docs/source_docs/tutorials/incompressible-fluid/steady-flow-around-sphere.rst
@@ -0,0 +1,248 @@
+.. _TutorialSteadyFlowAroundSphere:
+
+Steady flow around a sphere
+===========================
+
+This tutorial sets up steady, incompressible flow past an isolated sphere at
+:math:`Re=100`, well within the regime where the wake is steady and
+axisymmetric. It is the first tutorial in this series to use a fully
+three-dimensional embedded-boundary geometry -- both the
+:ref:`steady flow past a cylinder <TutorialSteadyFlowAroundCylinder>` and
+:ref:`backward-facing step <TutorialSteadyBackwardFacingStep>` tutorials are
+effectively two-dimensional, using a periodic or resolved-but-extruded
+spanwise direction. The sphere drag coefficient and wake recirculation
+length are computed as part of the analysis.
+
+Features
+--------
+
+* Incompressible fluid with no particles.
+* Steady-state advance.
+* Fully three-dimensional embedded-boundary sphere geometry.
+* Uniform freestream mass inflow, pressure outlet, and free-slip lateral
+  boundaries approximating an unconfined far field.
+* Plotfile output for velocity and pressure.
+* Surface-force (drag) reporting on the embedded boundary.
+* Comparison against the Schiller-Naumann drag correlation and the
+  Johnson & Patel (1999) :cite:p:`johnson99` / Taneda (1956)
+  :cite:p:`taneda56` recirculation-length data.
+* A four-level grid-refinement study.
+
+Case description
+-----------------
+
+The sphere has diameter :math:`D=1` (radius :math:`r=0.5`) and is centered
+at the origin. The domain extends 8 diameters upstream, 16 diameters
+downstream, and 8 diameters radially in :math:`y` and :math:`z`:
+
+.. math::
+
+   -8 \le x \le 16, \qquad -8 \le y \le 8, \qquad -8 \le z \le 8,
+
+giving a blockage ratio (frontal sphere area to domain cross-section) of
+about 0.3%, small enough that the lateral boundaries have a negligible
+effect on the resolved flow.
+
+The tutorial uses
+
+.. math::
+
+   \rho = 1,\qquad U_\infty = 1,\qquad D = 1,\qquad \mu = 0.01,
+
+giving
+
+.. math::
+
+   Re = \frac{\rho\, U_\infty\, D}{\mu} = 100.
+
+At this Reynolds number the wake is steady and axisymmetric about the flow
+direction; the transition to a steady but non-axisymmetric wake does not
+occur until :math:`Re \approx 212`, and vortex shedding does not begin until
+:math:`Re \approx 270`-:math:`300`.
+
+Important input sections
+--------------------------
+
+The sphere is defined through MFIX-Exa's built-in sphere embedded-boundary
+geometry:
+
+.. code-block:: text
+
+   mfix.geometry = "sphere"
+
+   sphere.radius = 0.5
+   sphere.center = 0.0  0.0  0.0
+   sphere.internal_flow = false
+
+Unlike the channel-flow cases in this series, the domain here represents an
+unconfined free stream rather than a physical enclosure, so the lateral
+boundaries use a free-slip condition instead of no-slip walls. The low-
+:math:`x` face is a uniform mass inflow, and the high-:math:`x` face is a
+pressure outlet:
+
+.. code-block:: text
+
+   bc.regions = inflow outflow side_ylo side_yhi side_zlo side_zhi
+
+   bc.inflow = mi
+   bc.inflow.fluid.density = 1.0
+   bc.inflow.fluid.inflow_type = velocity
+   bc.inflow.fluid.velocity = 1.0
+
+   bc.outflow = po
+   bc.outflow.pressure = 0.0
+
+   bc.side_ylo = slip
+   bc.side_yhi = slip
+   bc.side_zlo = slip
+   bc.side_zhi = slip
+
+No boundary condition entry is needed for the sphere itself -- embedded
+boundaries are no-slip by default. A separate region, slightly larger than
+the sphere, is used only to define where the drag report integrates:
+
+.. code-block:: text
+
+   regions.eb_sphere.shape = sphere
+   regions.eb_sphere.sphere.radius = 0.75
+   regions.eb_sphere.sphere.center = 0.0  0.0  0.0
+
+   mfix.reports.eb_drag.regions = eb_sphere
+   mfix.reports.eb_drag.eb_sphere.int = 500
+
+The fluid viscosity sets the Reynolds number:
+
+.. code-block:: text
+
+   fluid.viscosity.molecular.model = constant
+   fluid.viscosity.molecular.constant = 0.01
+
+As with the backward-facing-step case, the run uses the pseudo steady-state
+time advance:
+
+.. code-block:: text
+
+   mfix.steady_state = 1
+
+   mfix.steady_state_tol = 1.0e-6
+   mfix.steady_state_maxiter = 100000
+
+.. note::
+
+   The complete input file is located in the MFIX-Exa source directory as
+   ``tutorials/incompressible-fluid/inputs.steady-flow-past-sphere``
+
+Post-processing
+-----------------
+
+The run writes an AMReX plotfile named ``plt*`` every 5000 iterations into a
+directory, ``io``, in the run directory. Each plotfile saves the fluid
+velocity field, ``vel_g``, pressure, ``p_g``, and volume fraction,
+``volfrac``. The ``eb_drag`` report writes a time history of the surface
+force on the sphere, from which the drag coefficient
+
+.. math::
+
+   C_d = \frac{F_x}{\tfrac{1}{2}\rho U_\infty^2 A}, \qquad A = \pi r^2,
+
+is computed. All results were analyzed using a Jupyter notebook.
+
+A mid-plane slice through the wake, with the in-plane velocity streamlines
+overlaid, shows the expected steady, axisymmetric single recirculation
+bubble immediately behind the sphere:
+
+.. _fig_steady_flow_past_sphere_midplane:
+
+.. figure:: ./images/steady-flow-past-sphere-midplane_contour.png
+   :width: 80%
+   :align: center
+   :alt: Streamwise velocity and recirculation bubble behind a sphere at Reynolds number 100.
+
+   Streamwise velocity and recirculation bubble, :math:`Re=100`.
+
+Grid convergence
+------------------
+
+The case was run on four uniformly refined meshes, labeled by the number of
+cells in each direction:
+
+.. list-table:: Mesh resolutions used in the grid-refinement study
+   :header-rows: 1
+
+   * - Mesh
+     - Cells / diameter
+     - :math:`\Delta x = \Delta y = \Delta z`
+   * - 192x128x128
+     - 8
+     - 0.125000
+   * - 384x256x256
+     - 16
+     - 0.062500
+   * - 768x512x512
+     - 32
+     - 0.031250
+   * - 1536x1024x1024
+     - 64
+     - 0.015625
+
+Each finer level was initialized from the previous, coarser level's
+converged solution rather than started from a quiescent field, to reduce the
+iteration count needed to reach steady state.
+
+Two quantities are compared against reference data: the drag coefficient
+:math:`C_d`, against the Schiller-Naumann correlation
+(:math:`C_d = 1.0917` at :math:`Re=100`), and the recirculation length
+:math:`L/D`, measured from the rear stagnation point to the wake's zero-
+velocity crossing along the centerline, against Johnson & Patel
+:cite:p:`johnson99` and Taneda :cite:p:`taneda56`
+(:math:`L/D \approx 0.88`-:math:`0.89` at :math:`Re=100`).
+
+.. list-table:: Drag coefficient and recirculation length vs. mesh resolution
+   :header-rows: 1
+
+   * - Mesh
+     - :math:`C_d`
+     - % vs. Schiller-Naumann
+     - :math:`L/D`
+     - % vs. Johnson & Patel
+   * - 192x128x128
+     - 1.0469
+     - -4.10%
+     - 0.958
+     - +8.82%
+   * - 384x256x256
+     - 1.0396
+     - -4.78%
+     - 0.897
+     - +1.90%
+   * - 768x512x512
+     - 1.0299
+     - -5.66%
+     - 0.866
+     - -1.56%
+   * - 1536x1024x1024
+     - 1.0215
+     - -6.43%
+     - 0.868
+     - -1.32%
+   * - Reference
+     - 1.0917
+     -
+     - 0.88-0.89
+     -
+
+The recirculation length converges toward the Johnson & Patel / Taneda
+range over the first three mesh levels before ticking back up slightly at
+the finest resolution. The drag coefficient moves monotonically away from
+the Schiller-Naumann correlation with each refinement. Richardson
+extrapolation was attempted for both quantities across the available
+mesh-level triplets and did not produce a consistent order-of-accuracy
+estimate for either one; rather than present an extrapolated value that
+isn't actually trustworthy, the raw values above are reported as-is.
+
+
+.. note::
+
+   Because a stable grid-converged value could not be established, the
+   :math:`C_d` and :math:`L/D` figures above should be treated as
+   illustrative of the analysis workflow.
diff --git a/docs/source_docs/tutorials/index.rst b/docs/source_docs/tutorials/index.rst
index 50a096dd7bd59a8d14ba3de4e62ec845d585110f..2b616fe2c053adc637be63fdfb25e4a5315e53ee 100644
--- a/docs/source_docs/tutorials/index.rst
+++ b/docs/source_docs/tutorials/index.rst
@@ -7,4 +7,6 @@ used to build them.
 .. toctree::
    :maxdepth: 1
 
+   incompressible-fluid/steady-backward-facing-step
    incompressible-fluid/steady-flow-around-cylinder
+   incompressible-fluid/steady-flow-around-sphere