diff --git a/docs/developers_guide/ocean/api.md b/docs/developers_guide/ocean/api.md
index 185998b97..edf16d6c3 100644
--- a/docs/developers_guide/ocean/api.md
+++ b/docs/developers_guide/ocean/api.md
@@ -98,6 +98,7 @@
 
    init.Init
    init.Init.run
+   init.cosine_bell
 
    forward.Forward
    forward.Forward.compute_cell_count
diff --git a/docs/users_guide/ocean/tasks/cosine_bell.md b/docs/users_guide/ocean/tasks/cosine_bell.md
index b344fae14..18b1898d7 100644
--- a/docs/users_guide/ocean/tasks/cosine_bell.md
+++ b/docs/users_guide/ocean/tasks/cosine_bell.md
@@ -10,7 +10,7 @@ Bell test case as first described in
 but using the variant from Sec. 3a of
 [Skamarock and Gassmann](https://doi.org/10.1175/MWR-D-10-05056.1).  A flow
 field representing solid-body rotation transports a bell-shaped perturbation
-in a tracer $psi$ once around the sphere, returning to its initial location.
+in a tracer $\psi$ once around the sphere, returning to its initial location.
 
 The task is a convergence test with time step varying proportionately to grid
 size. The result of the `analysis` step of the task is a plot like the
@@ -133,6 +133,11 @@ sphere. `psi_0` and `radius`, $R$, are given as config options and may be
 altered by the user. In the `init` step we assign `debug_tracers_1`
 to $\psi$.
 
+```{image} images/cosine_bell_init.png
+:align: center
+:width: 500 px
+```
+
 The initial velocity is equatorial:
 
 $$
@@ -172,8 +177,8 @@ time_integrator = RK4
 rk4_dt_per_km = 3.0
 ```
 
-The convergence_eval_time, run duration and output interval are the period for 
-advection to make a full rotation around the globe, 24 days:
+The `convergence_eval_time`, `run_duration` and `output_interval` are the
+period for advection to make a full rotation around the globe, 24 days:
 
 ```cfg
 # config options for spherical convergence tests
diff --git a/docs/users_guide/ocean/tasks/images/cosine_bell_init.png b/docs/users_guide/ocean/tasks/images/cosine_bell_init.png
new file mode 100644
index 000000000..5a56a4252
Binary files /dev/null and b/docs/users_guide/ocean/tasks/images/cosine_bell_init.png differ
diff --git a/polaris/ocean/tasks/cosine_bell/init.py b/polaris/ocean/tasks/cosine_bell/init.py
index f1f316fda..251918a04 100644
--- a/polaris/ocean/tasks/cosine_bell/init.py
+++ b/polaris/ocean/tasks/cosine_bell/init.py
@@ -111,4 +111,18 @@ def run(self):
 
 
 def cosine_bell(max_value, ri, r):
+    """
+    Compute values according to cosine bell function
+
+    Parameters
+    ----------
+    max_value : float
+        Maximum value of the cosine bell function
+
+    ri : np.ndarray of type float
+        Distance from the center of the cosine bell in meters
+
+    r : float
+        Radius of the cosine bell in meters
+    """
     return max_value / 2.0 * (1.0 + np.cos(np.pi * np.divide(ri, r)))