diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 1963073423..4e5eb11ce5 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -11,9 +11,9 @@ version = "0.0.1"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "8bc0aaec0ca548eb6cf5f0d7d16351650c1ee956"
+git-tree-sha1 = "cad4c758c0038eea30394b1b671526921ca85b21"
 uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-version = "1.3.2"
+version = "1.4.0"
 weakdeps = ["ChainRulesCore"]
 
     [deps.AbstractFFTs.extensions]
@@ -142,9 +142,9 @@ version = "1.16.0"
 
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
-git-tree-sha1 = "9c209fb7536406834aa938fb149964b985de6c83"
+git-tree-sha1 = "02aa26a4cf76381be7f66e020a3eddeb27b0a092"
 uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
-version = "0.7.1"
+version = "0.7.2"
 
 [[deps.ColorBrewer]]
 deps = ["Colors", "JSON", "Test"]
@@ -205,9 +205,9 @@ version = "0.5.1"
 
 [[deps.ConcurrentUtilities]]
 deps = ["Serialization", "Sockets"]
-git-tree-sha1 = "96d823b94ba8d187a6d8f0826e731195a74b90e9"
+git-tree-sha1 = "5372dbbf8f0bdb8c700db5367132925c0771ef7e"
 uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb"
-version = "2.2.0"
+version = "2.2.1"
 
 [[deps.ConstructionBase]]
 deps = ["LinearAlgebra"]
@@ -276,10 +276,10 @@ uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 version = "0.27.25"
 
 [[deps.DocumenterCitations]]
-deps = ["Bibliography", "DataStructures", "Documenter", "Markdown", "Unicode"]
-git-tree-sha1 = "256f098ebde08780f058a9c84b4d08601b4ed2b4"
+deps = ["Bibliography", "Documenter", "Markdown", "OrderedCollections", "Unicode"]
+git-tree-sha1 = "b4b45d64eea45d96b5d35ec9ef4a0de15b7f2cdb"
 uuid = "daee34ce-89f3-4625-b898-19384cb65244"
-version = "0.2.12"
+version = "1.0.0"
 
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
@@ -355,9 +355,9 @@ uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
 [[deps.FillArrays]]
 deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
-git-tree-sha1 = "2250347838b28a108d1967663cba57bfb3c02a58"
+git-tree-sha1 = "e5556303fd8c9ad4a8fceccd406ef3433ddb4c45"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.3.0"
+version = "1.4.0"
 
 [[deps.FixedPointNumbers]]
 deps = ["Statistics"]
@@ -466,9 +466,9 @@ version = "1.12.2+2"
 
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
-git-tree-sha1 = "7f5ef966a02a8fdf3df2ca03108a88447cb3c6f0"
+git-tree-sha1 = "cb56ccdd481c0dd7f975ad2b3b62d9eda088f7e2"
 uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
-version = "1.9.8"
+version = "1.9.14"
 
 [[deps.HarfBuzz_jll]]
 deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
@@ -478,9 +478,9 @@ version = "2.8.1+1"
 
 [[deps.HypergeometricFunctions]]
 deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"]
-git-tree-sha1 = "0ec02c648befc2f94156eaef13b0f38106212f3f"
+git-tree-sha1 = "83e95aaab9dc184a6dcd9c4c52aa0dc26cd14a1d"
 uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
-version = "0.3.17"
+version = "0.3.21"
 
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
@@ -585,9 +585,9 @@ version = "1.0.0"
 
 [[deps.JLD2]]
 deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"]
-git-tree-sha1 = "42c17b18ced77ff0be65957a591d34f4ed57c631"
+git-tree-sha1 = "5df8278ad24772c0c6dbbeb97b162ccf29ced2a9"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-version = "0.4.31"
+version = "0.4.32"
 
 [[deps.JLLWrappers]]
 deps = ["Preferences"]
@@ -1233,18 +1233,18 @@ version = "0.1.1"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "Random", "StaticArraysCore"]
-git-tree-sha1 = "0da7e6b70d1bb40b1ace3b576da9ea2992f76318"
+git-tree-sha1 = "fffc14c695c17bfdbfa92a2a01836cdc542a1e46"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.6.0"
+version = "1.6.1"
 weakdeps = ["Statistics"]
 
     [deps.StaticArrays.extensions]
     StaticArraysStatisticsExt = "Statistics"
 
 [[deps.StaticArraysCore]]
-git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
+git-tree-sha1 = "1d5708d926c76a505052d0d24a846d5da08bc3a4"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.0"
+version = "1.4.1"
 
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
@@ -1279,9 +1279,9 @@ version = "1.3.0"
 
 [[deps.StringEncodings]]
 deps = ["Libiconv_jll"]
-git-tree-sha1 = "33c0da881af3248dafefb939a21694b97cfece76"
+git-tree-sha1 = "b765e46ba27ecf6b44faf70df40c57aa3a547dcb"
 uuid = "69024149-9ee7-55f6-a4c4-859efe599b68"
-version = "0.3.6"
+version = "0.3.7"
 
 [[deps.StructArrays]]
 deps = ["Adapt", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 601001d60e..01b3908b26 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,4 +12,5 @@ TimesDates = "bdfc003b-8df8-5c39-adcd-3a9087f5df4a"
 
 [compat]
 CairoMakie = "0.10 - 0.10.4" # v0.10.5, v0.10.6 very slow with contourf + prescribed levels
+DocumenterCitations = "1"
 Literate = "≥2.9.0"
diff --git a/docs/make.jl b/docs/make.jl
index 0a53ee49ef..44b793894d 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -31,7 +31,7 @@ end
     using Oceananigans.BoundaryConditions: Flux, Value, Gradient, Open
 
     bib_filepath = joinpath(dirname(@__FILE__), "oceananigans.bib")
-    bib = CitationBibliography(bib_filepath)
+    bib = CitationBibliography(bib_filepath, style=:authoryear)
 
     #####
     ##### Generate examples
@@ -170,7 +170,8 @@ pages = [
 format = Documenter.HTML(collapselevel = 1,
                          prettyurls = get(ENV, "CI", nothing) == "true",
                          canonical = "https://clima.github.io/OceananigansDocumentation/stable/",
-                         mathengine = MathJax3())
+                         mathengine = MathJax3(),
+                         assets = String["assets/citations.css"])
 
 makedocs(bib, sitename = "Oceananigans.jl",
               authors = "Climate Modeling Alliance and contributors",
diff --git a/docs/oceananigans.bib b/docs/oceananigans.bib
index 35f7ff6802..5111f6d167 100644
--- a/docs/oceananigans.bib
+++ b/docs/oceananigans.bib
@@ -88,7 +88,7 @@ @book{Hesthaven07
 }
 
 @article{Marshall97FV,
-	title = {A finite-volume, incompressible {Navier--Stokes} model for studies of the ocean on parallel computers},
+	title = {A finite-volume, incompressible Navier--Stokes model for studies of the ocean on parallel computers},
 	volume = {102},
 	doi = {10.1029/96JC02775},
 	number = {C3},
@@ -176,9 +176,9 @@ @book{Boussinesq1903
 }
 
 @book{Boussinesq1877,
-  title={Essai sur la th{\'e}orie des eaux courantes},
+  title={Essai sur la théorie des eaux courantes},
   author={Boussinesq, J.},
-  series={M{\'e}moires pr{\'e}sent{\'e}s par divers savants {\`a} l'Acad{\'e}mie des sciences de l'Institut national de France},
+  series={Mémoires présentés par divers savants à l'Académie des sciences de l'Institut national de France},
   year={1877},
   publisher={Impr. Nationale}
 }
@@ -219,7 +219,7 @@ @article{Orszag86
 }
 
 @article{Brown01,
-	title = {Accurate projection methods for the incompressible {Navier--Stokes} equations},
+	title = {Accurate projection methods for the incompressible Navier--Stokes equations},
 	volume = {168},
 	doi = {10/djz53c},
 	number = {2},
@@ -253,7 +253,7 @@ @article{Harlow65
 
 
 @book{Kumar16,
-  title = {Local {BVP} methods for the computation of cell-face velocities in the incompressible {Navier--Stokes} equations},
+  title = {Local {BVP} methods for the computation of cell-face velocities in the incompressible Navier--Stokes equations},
   author = {N. Kumar and Boonkkamp, J.H.M. and B. Koren},
   year = {2016},
   series = {CASA-report},
@@ -262,7 +262,7 @@ @book{Kumar16
 }
 
 @article{Sani81,
-	title = {The cause and cure (?) of the spurious pressures generated by certain {FEM} solutions of the incompressible {Navier--Stokes} equations: {Part} 1},
+	title = {The cause and cure (?) of the spurious pressures generated by certain {FEM} solutions of the incompressible Navier--Stokes equations: {Part} 1},
 	volume = {1},
 	doi = {10.1002/fld.1650010104},
 	number = {1},
@@ -643,7 +643,7 @@ @article{Verstappen18
 	title = {How much eddy dissipation is needed to counterbalance the nonlinear production of small, unresolved scales in a large-eddy simulation of turbulence?},
 	volume = {176},
 	doi = {10.1016/j.compfluid.2016.12.016},
-	journal = {Computers \& Fluids},
+	journal = {Computers & Fluids},
 	author = {Verstappen, Roel},
 	year = {2018},
 	pages = {276--284}
@@ -664,7 +664,7 @@ @article{Bruneau06
 	volume = {35},
 	doi = {10/fvshbk},
 	number = {3},
-	journal = {Computers \& Fluids},
+	journal = {Computers & Fluids},
 	author = {Bruneau, Charles-Henri and Saad, Mazen},
 	year = {2006},
 	pages = {326--348}
@@ -675,7 +675,7 @@ @article{Botella98
 	volume = {27},
 	doi = {10/chbd5b},
 	number = {4},
-	journal = {Computers \& Fluids},
+	journal = {Computers & Fluids},
 	author = {Botella, O. and Peyret, R.},
 	year = {1998},
 	pages = {421--433}
diff --git a/docs/src/appendix/fractional_step.md b/docs/src/appendix/fractional_step.md
index 58cf3faf44..90f396fbff 100644
--- a/docs/src/appendix/fractional_step.md
+++ b/docs/src/appendix/fractional_step.md
@@ -30,7 +30,7 @@ problem for the pressure ``p^{n+1}`` with the boundary condition
   \boldsymbol{\hat{n}} \boldsymbol{\cdot} \boldsymbol{\nabla} p^{n+1} |_{\partial\Omega} = 0 \, .
 ```
 
-[Orszag86](@cite) and [Brown01](@cite) raise an important issue regarding these fractional step 
+[Orszag86](@citet) and [Brown01](@citet) raise an important issue regarding these fractional step 
 methods, which is that "*while the velocity can be reliably computed to second-order accuracy 
 in time and space, the pressure is typically only first-order accurate in the ``L_\infty``-norm.*" 
 The numerical boundary conditions must be carefully accounted for to ensure the second-order 
diff --git a/docs/src/appendix/staggered_grid.md b/docs/src/appendix/staggered_grid.md
index 0a594a1562..0b4da37aa5 100644
--- a/docs/src/appendix/staggered_grid.md
+++ b/docs/src/appendix/staggered_grid.md
@@ -18,7 +18,7 @@ on the top and bottom edges of the pressure control volumes. The indexing conven
 
 This staggered arrangement of variables is more complicated than the collocated grid arrangement but is greatly
 beneficial as it avoids the odd-even decoupling between the pressure and velocity if they are stored at the same
-positions. §6.1 of [Patankar80](@cite) discusses this problem in the presence of a zigzag pressure field: on a 1D
+positions. §6.1 of [Patankar80](@citet) discusses this problem in the presence of a zigzag pressure field: on a 1D
 collocated grid the velocity at the point ``i`` is influenced by the pressure at points ``i-1`` and ``i+1``, and a zigzag
 pressure field will be felt as a uniform pressure, which is obviously wrong and would reduce the accuracy of the
 solution. The pressure is effectively taken from a coarser grid than what is actually used. The basic problem is that
@@ -29,23 +29,23 @@ From the viewpoint of linear algebra, these spurious pressure modes correspond t
 pressure projection operator with eigenvalue zero and are thus indistinguishable from a uniform pressure field
 [Sani81](@cite).
 
-The staggered grid was first introduced by [Harlow65](@cite) with their *marker and cell* method. In meteorology
-and oceanography, this particular staggered grid configuration is referred to as the Arakawa C-grid after [Arakawa77](@cite), who
+The staggered grid was first introduced by [Harlow65](@citet) with their *marker and cell* method. In meteorology
+and oceanography, this particular staggered grid configuration is referred to as the Arakawa C-grid after [Arakawa77](@citet), who
 investigated four different staggered grids and the unstaggered A-grid for use in an atmospheric model.
 
-[Arakawa77](@cite) investigated the dispersion relation of inertia-gravity waves[^2] traveling in the ``x``-direction
+[Arakawa77](@citet) investigated the dispersion relation of inertia-gravity waves[^2] traveling in the ``x``-direction
 ```math
   \omega^2 = f^2 + gHk^2 \, ,
 ```
 in the linearized rotating shallow-water equations for five grids. Here ``\omega`` is the angular frequency, ``H`` is the
 height of the fluid and ``k`` is the wavenumber in the ``x``-direction. Looking at the effect of spatial discretization
 error on the frequency of these waves they find that the B and C-grids reproduce the dispersion relation most closely
-out of the five [Arakawa77](@cite) (Figure 5). In particular, the dispersion relation for the C-grid is given by
+out of the five [Arakawa77](@citet) (Figure 5). In particular, the dispersion relation for the C-grid is given by
 ```math
   \omega^2 = f^2 \left[ \cos^2 \left( \frac{k\Delta}{2} \right)
              + 4 \left( \frac{\lambda}{\Delta} \right)^2 \sin^2 \left( \frac{k\Delta}{2} \right) \right] \, ,
 ```
-where ``\lambda`` is the wavelength and ``\Delta`` is the grid spacing. Paraphrasing p. 184 of [Arakawa77](@cite): The
+where ``\lambda`` is the wavelength and ``\Delta`` is the grid spacing. Paraphrasing p. 184 of [Arakawa77](@citet): The
 wavelength of the shortest resolvable wave is ``2\Delta`` with corresponding wavenumber ``k = \pi/\Delta`` so it is
 sufficient to evaluate the dispersion relation over the range ``0 < k \Delta < \pi``. The frequency is monotonically
 increasing for ``\lambda / \Delta > \frac{1}{2}`` and monotonically decreasing for ``\lambda / \Delta < \frac{1}{2}``. For the
@@ -55,8 +55,8 @@ oscillations or stationary waves, which is bad.
 
 The B and C-grids are less oscillatory than the others and quite faithfully simulate geostrophic adjustment. However,
 the C-grid is the only one that faithfully reproduces the two-dimensional dispersion relation ``\omega^2(k, \ell)``, all
-the other grids have false maxima, and so [Arakawa77](@cite) conclude that the C-grid is best for simulating geostrophic
+the other grids have false maxima, and so [Arakawa77](@citet) conclude that the C-grid is best for simulating geostrophic
 adjustment except for abnormal situations in which ``\lambda / \Delta`` is less than or close to 1. This seems to have held
 true for most atmospheric and oceanographic simulations as the C-grid is popular and widely used.
 
-[^2]: Apparently also called Poincaré waves, Sverdrup waves, and *rotational gravity waves* §13.9 of [Kundu15](@cite).
+[^2]: Apparently also called Poincaré waves, Sverdrup waves, and *rotational gravity waves* §13.9 of [Kundu15](@citet).
diff --git a/docs/src/assets/citations.css b/docs/src/assets/citations.css
new file mode 100644
index 0000000000..844fd5ec85
--- /dev/null
+++ b/docs/src/assets/citations.css
@@ -0,0 +1,18 @@
+.citation dl {
+    display: grid;
+    grid-template-columns: max-content auto; }
+  .citation dt {
+    grid-column-start: 1; }
+  .citation dd {
+    grid-column-start: 2;
+    margin-bottom: 0.75em; }
+  .citation ul {
+   padding: 0 0 2.25em 0;
+   margin: 0;
+   list-style: none;}
+  .citation ul li {
+   text-indent: -2.25em;
+   margin: 0.33em 0.5em 0.5em 2.25em;}
+  .citation ol li {
+   padding-left:0.75em;}
+  
\ No newline at end of file
diff --git a/docs/src/model_setup/buoyancy_and_equation_of_state.md b/docs/src/model_setup/buoyancy_and_equation_of_state.md
index c7aad171b8..737a98965d 100644
--- a/docs/src/model_setup/buoyancy_and_equation_of_state.md
+++ b/docs/src/model_setup/buoyancy_and_equation_of_state.md
@@ -185,7 +185,7 @@ SeawaterBuoyancy{Float64}:
 ### Idealized nonlinear equations of state
 
 Instead of a linear equation of state, six idealized (second-order) nonlinear equations of state
-as described by [Roquet15Idealized](@cite) may be used. These equations of state are provided
+as described by [Roquet15Idealized](@citet) may be used. These equations of state are provided
 via the [SeawaterPolynomials.jl](https://github.com/CliMA/SeawaterPolynomials.jl) package.
 
 ```jldoctest buoyancy
@@ -208,7 +208,7 @@ SeawaterBuoyancy{Float64}:
 ### TEOS-10 equation of state
 
 A high-accuracy 55-term polynomial approximation to the TEOS-10 equation of state suitable for use in
-Boussinesq models as described by [Roquet15TEOS](@cite) is implemented in the
+Boussinesq models as described by [Roquet15TEOS](@citet) is implemented in the
 [SeawaterPolynomials.jl](https://github.com/CliMA/SeawaterPolynomials.jl) package and may be used.
 
 ```jldoctest buoyancy
diff --git a/docs/src/numerical_implementation/elliptic_solvers.md b/docs/src/numerical_implementation/elliptic_solvers.md
index 0c483685f5..7ade5265f2 100644
--- a/docs/src/numerical_implementation/elliptic_solvers.md
+++ b/docs/src/numerical_implementation/elliptic_solvers.md
@@ -27,8 +27,8 @@ Discretizing elliptic problems that can be solved via a classical separation-of-
 equation, results in a linear system of equations ``M \boldsymbol{x} = \boldsymbol{y}`` where ``M`` is a real symmetric matrix of block
 tridiagonal form. This allows for the matrix to be decomposed and solved efficiently, provided that the eigenvalues and
 eigenvectors of the blocks are known (§2) [Buzbee70](@cite). In the case of Poisson's equation on a rectangle,
-[Hockney65](@cite) has taken advantage of the fact that the fast Fourier transform can be used to perform the matrix
-multiplication steps resulting in an even more efficient method. [Schumann88](@cite) describe the implementation of such
+[Hockney65](@citet) has taken advantage of the fact that the fast Fourier transform can be used to perform the matrix
+multiplication steps resulting in an even more efficient method. [Schumann88](@citet) describe the implementation of such
 an algorithm for Poisson's equation on a staggered grid with Dirichlet, Neumann, and periodic boundary conditions.
 
 The method can be explained easily by taking the Fourier transform of both sides of \eqref{eq:poisson-pressure} to yield
@@ -51,7 +51,7 @@ coefficients of the solution are easily solved for
 \widehat{p}_{NH}(i, j, k) = - \frac{\widehat{\mathscr{F}}(i, j, k)}{\lambda^x_i + \lambda^y_j + \lambda^z_k} \, .
 ```
 
-The eigenvalues are given by [Schumann88](@cite) and can also be tediously derived by plugging in the definition of the
+The eigenvalues are given by [Schumann88](@citet) and can also be tediously derived by plugging in the definition of the
 discrete Fourier transform into \eqref{eq:poisson-spectral}:
 ```math
 \begin{align}
@@ -81,8 +81,8 @@ system along one of the dimensions and utilizing cyclic reduction. This results
 reduction* or ``\text{FACR}(\ell)`` algorithm (with ``\ell`` cyclic reduction steps) which requires only
 ``\mathcal{O}(N \log_2\log_2 N)`` operations provided the optimal number of cyclic reduction steps is taken, which is
 ``\ell = \log_2 \log_2 n`` where ``n`` is the number of grid points in the cyclic reduction dimension. The FACR algorithm
-was first developed by [Hockney69](@cite) and is well reviewed by [Swarztrauber77](@cite) then further benchmarked and
-extended by [Temperton79](@cite) and [Temperton80](@cite).
+was first developed by [Hockney69](@citet) and is well reviewed by [Swarztrauber77](@citet) then further benchmarked and
+extended by [Temperton79](@citet) and [Temperton80](@citet).
 
 Furthermore, the FACR algorithm removes the restriction that the grid is uniform in one of the dimensions so it can
 be utilized to implement a fast Poisson solver for vertically stretched grids if the cyclic reduction is applied in the
@@ -113,7 +113,7 @@ Discretizing the ``\partial_z^2`` derivative and equating the term inside the br
 ## Cosine transforms on the GPU
 
 Unfortunately cuFFT does not provide cosine transforms and so we must write our own fast cosine 
-transforms for the GPU. We implemented the fast 1D and 2D cosine transforms described by [Makhoul80](@cite) 
+transforms for the GPU. We implemented the fast 1D and 2D cosine transforms described by [Makhoul80](@citet) 
 which compute it by applying the regular Fourier transform to a permuted version of the array.
 
 In this section we will be using the DCT-II as the definition of the forward cosine transform 
diff --git a/docs/src/numerical_implementation/large_eddy_simulation.md b/docs/src/numerical_implementation/large_eddy_simulation.md
index 8f6597ff48..a19504b3e5 100644
--- a/docs/src/numerical_implementation/large_eddy_simulation.md
+++ b/docs/src/numerical_implementation/large_eddy_simulation.md
@@ -5,25 +5,25 @@ sub-grid scale motions. This is done usually be assuming eddy viscosity and eddy
 estimate for the eddy viscosity ``\nu_e`` and diffusivity ``\kappa_e``.
 
 Much of the early work on LES was motivated by the study of atmospheric boundary layer turbulence, being developed
-by [Smagorinsky63](@cite) and [Lilly66](@cite), then first implemented by [Deardorff70](@cite) and [Deardorff74](@cite).
+by [Smagorinsky63](@citet) and [Lilly66](@citet), then first implemented by [Deardorff70](@citet) and [Deardorff74](@citet).
 
-In the LES framework, the Navier-Stokes equations are averaged in the same way as [Reynolds1895](@cite) except that the
+In the LES framework, the Navier-Stokes equations are averaged in the same way as [Reynolds1895](@citet) except that the
 mean field ``\overline{\boldsymbol{v}}`` is obtained via convolution with a filter convolution kernel ``G``
 ```math
 \overline{\boldsymbol{v}(\boldsymbol{x}, t)} = G \star \boldsymbol{v} =
   \int_{-\infty}^\infty \int_{-\infty}^\infty
   \boldsymbol{v}(\boldsymbol{x}^\prime, t) G(\boldsymbol{x} - \boldsymbol{x}^\prime, t - \tau) \, \mathrm{d}\boldsymbol{x}^\prime \, \mathrm{d} \tau \, ,
 ```
-as described by [Leonard75](@cite) who introduced the general filtering formalism.
+as described by [Leonard75](@citet) who introduced the general filtering formalism.
 
 The ``\overline{v_i^\prime v_j^\prime}`` terms are now components of what is called the sub-grid scale (SGS) stress
 tensor ``\tau^\text{SGS}_{ij}``, which looks the same as the Reynolds stress tensor so we will drop the SGS superscript.
 
 It is probably important to note that the large eddy simulation filtering operation does not satisfy the properties
-of a Reynolds operator (§2.1)[sagaut06](@cite) and that in general, the filtered residual is not zero:
+of a Reynolds operator (§2.1) [sagaut06](@cite) and that in general, the filtered residual is not zero:
 ``\overline{\boldsymbol{v}^\prime(\boldsymbol{x}, t)} \ne 0``.
 
-§13.2 of [Pope00](@cite) lists a number of popular choices for the filter function ``G``. For practical reasons we
+§13.2 of [Pope00](@citet) lists a number of popular choices for the filter function ``G``. For practical reasons we
 simply employ the box kernel
 ```math
   \begin{equation}
@@ -32,7 +32,7 @@ simply employ the box kernel
   \end{equation}
 ```
 where ``H`` is the Heaviside function, ``\Delta`` is the grid spacing, and ``t_n`` is the current time step. With
-\eqref{eq:box-kernel} we get back the averaging operator originally used by [Deardorff70](@cite)
+\eqref{eq:box-kernel} we get back the averaging operator originally used by [Deardorff70](@citet)
 ```math
 \overline{\boldsymbol{v}(x, y, z, t)} =
   \frac{1}{\Delta x \Delta y \Delta z}
@@ -46,7 +46,7 @@ which if evaluated at the cell centers just returns the cell averages we already
 
 ## Smagorinsky-Lilly model
 
-[Smagorinsky63](@cite) estimated the eddy viscosity ``\nu_e`` via a characteristic length scale ``\Delta`` times a velocity
+[Smagorinsky63](@citet) estimated the eddy viscosity ``\nu_e`` via a characteristic length scale ``\Delta`` times a velocity
 scale given by ``\Delta |\overline{S}|`` where ``|\overline{S}| = \sqrt{2\overline{S}_{ij}\overline{S}_{ij}}``. Thus the
 SGS stress tensor is given by
 ```math
@@ -57,7 +57,7 @@ The eddy diffusivities are calculated via ``\kappa_e = \nu_e / \text{Pr}_t`` whe
 ``\text{Pr}_t`` is usually chosen to be ``\mathcal{O}(1)`` from experimental observations.
 
 Assuming that the SGS energy cascade is equal to the overall dissipation rate ``\varepsilon`` from the
-[Kolmogorov41](@cite) theory, [Lilly66](@cite) was able to derive a value of
+[Kolmogorov41](@citet) theory, [Lilly66](@citet) was able to derive a value of
 ```math
 C_s = \left( \frac{3}{2}C_K\pi^\frac{4}{3} \right)^{-\frac{3}{4}} \approx 0.16 \, ,
 ```
@@ -71,16 +71,16 @@ convection.
 ## Anisotropic minimum dissipation models
 
 Minimum-dissipation eddy-viscosity models are a class of LES closures that use the minimum eddy dissipation required to
-dissipate the energy of sub-grid scale motion. [Rozema15](@cite) proposed the first minimum-dissipation model
+dissipate the energy of sub-grid scale motion. [Rozema15](@citet) proposed the first minimum-dissipation model
 appropriate for use on anisotropic grids, termed the *anisotropic minimum dissipation* (AMD) model.
 
 It has a number of desirable properties over Smagorinsky-type closures: it is more cost-effective than dynamic
 Smagorinsky, it appropriately switches off in laminar and transitional flows, and it is consistent with the exact SGS
-stress tensor on both isotropic and anisotropic grids. [Abkar16](@cite) extended the AMD model to model SGS scalar
-fluxes for tracer transport. [Abkar17](@cite) further extended the model to include a buoyancy term that accounts for
+stress tensor on both isotropic and anisotropic grids. [Abkar16](@citet) extended the AMD model to model SGS scalar
+fluxes for tracer transport. [Abkar17](@citet) further extended the model to include a buoyancy term that accounts for
 the contribution of buoyant forces to the production and suppression of turbulence.
 
-[Vreugdenhil18](@cite) derive a modified AMD model by following the requirement suggested by [Verstappen18](@cite),
+[Vreugdenhil18](@citet) derive a modified AMD model by following the requirement suggested by [Verstappen18](@citet),
 which entail normalising the displacement, the velocity, and the velocity gradient by the filter width to ensure that
 the resulting eddy dissipation properly counteracts the spurious kinetic energy transferred by convective nonlinearity,
 to derive a modified AMD model.
@@ -132,12 +132,12 @@ so that the normalized rate of strain tensor is
 ```
 
 In equations \eqref{eq:nu-dagger}--\eqref{eq:S-hat}, ``C`` is a modified Poincaré "constant" that is independent from
-the filter width ``\Delta`` but does depend on the accuracy of the discretization method used. [Abkar16](@cite) cite
+the filter width ``\Delta`` but does depend on the accuracy of the discretization method used. [Abkar16](@citet) cite
 ``C^2 = \frac{1}{12}`` for a spectral method and ``C^2 = \frac{1}{3}`` for a second-order accurate scheme. ``\Delta_i`` is
 the filter width in the ``x_i``-direction, and ``\Delta`` is given by the square root of the harmonic mean of the squares
 of the filter widths in each direction
 ```math
     \frac{1}{\Delta^2} = \frac{1}{3} \left( \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2} \right) \, .
 ```
-The term multiplying ``C_b`` is the buoyancy modification introduced by [Abkar17](@cite) and is small for weakly
+The term multiplying ``C_b`` is the buoyancy modification introduced by [Abkar17](@citet) and is small for weakly
 stratified flows. We have introduced the ``C_b`` constant so that the buoyancy modification term may be turned on and off.
diff --git a/docs/src/numerical_implementation/spatial_operators.md b/docs/src/numerical_implementation/spatial_operators.md
index 3878ff6a1d..64b5ba8618 100644
--- a/docs/src/numerical_implementation/spatial_operators.md
+++ b/docs/src/numerical_implementation/spatial_operators.md
@@ -2,7 +2,7 @@
 
 To calculate the various terms and perform the time-stepping, discrete difference and interpolation 
 operators must be designed from which all the terms, such as momentum advection and Laplacian 
-diffusion, may be constructed. Much of the material in this section is derived from [Marshall97FV](@cite).
+diffusion, may be constructed. Much of the material in this section is derived from [Marshall97FV](@citet).
 
 ## Differences
 
diff --git a/docs/src/numerical_implementation/time_stepping.md b/docs/src/numerical_implementation/time_stepping.md
index 218c42c473..6831eea1b2 100644
--- a/docs/src/numerical_implementation/time_stepping.md
+++ b/docs/src/numerical_implementation/time_stepping.md
@@ -70,7 +70,7 @@ for any time-dependent function ``G(t)``, while a second-order Adams-Bashforth m
         - \left ( \tfrac{1}{2} + \chi \right ) G^{n-1} \right ] \, ,
     \end{equation}
 ```
-where ``\chi`` is a parameter. [Ascher95](@cite) claim that ``\chi = \tfrac{1}{8}`` is optimal; 
+where ``\chi`` is a parameter. [Ascher95](@citet) claim that ``\chi = \tfrac{1}{8}`` is optimal; 
 ``\chi = -\tfrac{1}{2}`` yields the forward Euler scheme.
 
 Combining the equation \eqref{eq:intermediate-velocity-field} for ``\boldsymbol{v}^\star`` and the time integral
diff --git a/docs/src/numerical_implementation/turbulence_closures.md b/docs/src/numerical_implementation/turbulence_closures.md
index b9e4cf124b..693efd8124 100644
--- a/docs/src/numerical_implementation/turbulence_closures.md
+++ b/docs/src/numerical_implementation/turbulence_closures.md
@@ -1,11 +1,11 @@
 # [Turbulence closures](@id numerical_closures)
 
 To truly simulate and resolve turbulence at high Reynolds number (so basically all interesting flows) would require
-you resolve all motions down to the [Kolmogorov41](@cite) length scale ``\eta = (\nu^3 / \varepsilon)^{1/4}`` where
+you resolve all motions down to the [Kolmogorov41](@citet) length scale ``\eta = (\nu^3 / \varepsilon)^{1/4}`` where
 ``\nu`` is the kinematic viscosity and ``\varepsilon`` the average rate of dissipation of turbulence kinetic energy per
 unit mass.
 
-As pointed out way back by [Corrsin61](@cite), to run a simulation on a horizontal domain about 10 times the size of an
+As pointed out way back by [Corrsin61](@citet), to run a simulation on a horizontal domain about 10 times the size of an
 "average eddy" with 100 vertical levels and where the grid spacing is given by ``\eta`` would require the computer to
 store on the order of ``10^{14}`` variables.[^1] This is still impractical today, although may be within
 reach in less than a decade. He ends by suggesting the use of an analog rather digital computer---a tank of water.
@@ -22,9 +22,9 @@ To have any hope of simulating high Reynolds number flows we need some way of re
 
 ## Reynolds-averaged Navier–Stokes equations
 
-Following [Reynolds1895](@cite) we can decompose flow variables such as velocity ``\boldsymbol{v}`` into the mean component
+Following [Reynolds1895](@citet), we can decompose flow variables such as velocity ``\boldsymbol{v}`` into the mean component
 ``\overline{\boldsymbol{v}}`` and the fluctuating component ``\boldsymbol{v}^\prime`` so that ``\boldsymbol{v} = \overline{\boldsymbol{v}} + \boldsymbol{v}^\prime``
-[see §4 of [Pope00](@cite) for a modern discussion].
+[see §4 of [Pope00](@citet) for a modern discussion].
 
 Expressing the Navier-Stokes equations in tensor notation
 ```math
@@ -62,7 +62,7 @@ This is kind of hopeless so we will have to find some way to model the Reynolds
 
 ## Gradient-diffusion hypothesis and eddy viscosity models
 
-The *gradient-diffusion hypothesis*, due to [Boussinesq1877](@cite), assumes that the transport of scalar fluxes
+The *gradient-diffusion hypothesis*, due to [Boussinesq1877](@citet), assumes that the transport of scalar fluxes
 such as ``\overline{\boldsymbol{v}^\prime c^\prime}`` and ``\overline{v_i^\prime v_j^\prime}`` occurs down the mean scalar gradient
 ``\boldsymbol{\nabla} \overline{c}`` as if they are being diffused (§4.4) [Pope00](@cite). This is in analogy with how momentum transfer by
 molecular motion in a gas can be described by a molecular viscosity.
diff --git a/docs/src/physics/buoyancy_and_equations_of_state.md b/docs/src/physics/buoyancy_and_equations_of_state.md
index 6333b00b07..a695602412 100644
--- a/docs/src/physics/buoyancy_and_equations_of_state.md
+++ b/docs/src/physics/buoyancy_and_equations_of_state.md
@@ -33,4 +33,4 @@ and ``\beta`` is the haline contraction coefficient.
 
 ### Nonlinear equation of state
 
-Buoyancy is determined by the simplified equations of state introduced by [Roquet15TEOS](@cite).
+Buoyancy is determined by the simplified equations of state introduced by [Roquet15TEOS](@citet).
diff --git a/docs/src/physics/coriolis_forces.md b/docs/src/physics/coriolis_forces.md
index 77d58f74c1..e71eddffa7 100644
--- a/docs/src/physics/coriolis_forces.md
+++ b/docs/src/physics/coriolis_forces.md
@@ -59,5 +59,5 @@ the locally vertical and the locally horizontal components of the rotation vecto
     \boldsymbol{f} = \left[ 2\Omega\cos\varphi_0 \left( 1 -  \frac{z}{R} \right) + \gamma y \right] \boldsymbol{\hat y}
            + \left[ 2\Omega\sin\varphi_0 \left( 1 + 2\frac{z}{R} \right) + \beta  y \right] \boldsymbol{\hat z} \, ,
 ```
-as can be found in the paper by [DellarJFM2011](@cite) where 
+as can be found in the paper by [DellarJFM2011](@citet) where 
 ``\beta = 2 \Omega \cos \varphi_0 / R`` and ``\gamma = -4 \Omega \sin \varphi_0 / R``.
diff --git a/docs/src/physics/turbulence_closures.md b/docs/src/physics/turbulence_closures.md
index 0ab21dec9f..961bdbf463 100644
--- a/docs/src/physics/turbulence_closures.md
+++ b/docs/src/physics/turbulence_closures.md
@@ -65,7 +65,7 @@ diffusivity components ``\kappa_h`` and ``\kappa_z``.
 
 ## Smagorinsky-Lilly turbulence closure
 
-In the turbulence closure proposed by [Lilly62](@cite) and [Smagorinsky63](@cite),
+In the turbulence closure proposed by [Lilly62](@citet) and [Smagorinsky63](@citet),
 the subgrid stress associated with unresolved turbulent motions is modeled diffusively via
 ```math
 \tau_{ij} = - 2 \nu_e \Sigma_{ij} \, ,
@@ -110,8 +110,8 @@ Both ``Pr`` and ``\kappa`` may be set independently for each tracer.
 
 ## Anisotropic minimum dissipation (AMD) turbulence closure
 
-The anisotropic minimum dissipation (AMD) model proposed by [Verstappen18](@cite) and was
-described and tested by [Vreugdenhil18](@cite). The AMD model uses an eddy diffusivity hypothesis
+The anisotropic minimum dissipation (AMD) model proposed by [Verstappen18](@citet) and was
+described and tested by [Vreugdenhil18](@citet). The AMD model uses an eddy diffusivity hypothesis
 similar the Smagorinsky-Lilly model. In the AMD model, the eddy viscosity and diffusivity for each
 tracer are defined in terms of eddy viscosity and diffusivity *predictors* ``\nu_e^\dagger`` and ``\kappa_e^\dagger``,
 such that
diff --git a/paper/paper.bib b/paper/paper.bib
index 7796fddfa4..18a46a4a85 100644
--- a/paper/paper.bib
+++ b/paper/paper.bib
@@ -1,5 +1,5 @@
 @article{Marshall1997,
-  title   = {A finite-volume, incompressible {Navier} {Stokes} model for studies of the ocean on parallel computers},
+  title   = {A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers},
   author  = {Marshall, John and Adcroft, Alistair and Hill, Chris and Perelman, Lev and Heisey, Curt},
   journal = {Journal of Geophysical Research: Oceans},
   volume  = {102},
diff --git a/src/Solvers/discrete_transforms.jl b/src/Solvers/discrete_transforms.jl
index e6e06085b8..744ea60d3d 100644
--- a/src/Solvers/discrete_transforms.jl
+++ b/src/Solvers/discrete_transforms.jl
@@ -41,7 +41,7 @@ twiddle_factors(arch, grid, dim) = nothing
 """
     twiddle_factors(arch::GPU, grid, dims)
 
-Twiddle factors are needed to perform DCTs on the GPU. See equations (19a) and (22) of [Makhoul80](@cite)
+Twiddle factors are needed to perform DCTs on the GPU. See equations (19a) and (22) of [Makhoul80](@citet)
 for the forward and backward twiddle factors respectively.
 """
 function twiddle_factors(arch::GPU, grid, dims)
diff --git a/src/Solvers/index_permutations.jl b/src/Solvers/index_permutations.jl
index ce53a753af..a7599289c8 100644
--- a/src/Solvers/index_permutations.jl
+++ b/src/Solvers/index_permutations.jl
@@ -13,7 +13,7 @@ Permute `i` such that, for example, `i ∈ 1:N` becomes
 
 for `N=8` and `N=9` respectively.
 
-See equation (20) of [Makhoul80](@cite).
+See equation (20) of [Makhoul80](@citet).
 """
 @inline permute_index(i, N)::Int = ifelse(isodd(i),
                                           Base.unsafe_trunc(Int, i/2) + 1,
@@ -31,7 +31,7 @@ for example, `i ∈ 1:N` becomes
 
 for `N=8` and `N=9` respectively.
 
-See equation (20) of [Makhoul80](@cite).
+See equation (20) of [Makhoul80](@citet).
 """
 @inline unpermute_index(i, N) = ifelse(i <= ceil(N/2), 2i-1, 2(N-i+1))
 
diff --git a/src/TimeSteppers/runge_kutta_3.jl b/src/TimeSteppers/runge_kutta_3.jl
index 397f0d4167..1f61df1056 100644
--- a/src/TimeSteppers/runge_kutta_3.jl
+++ b/src/TimeSteppers/runge_kutta_3.jl
@@ -27,7 +27,7 @@ end
 Return a 3rd-order Runge0Kutta timestepper (`RungeKutta3TimeStepper`) on `grid` and with `tracers`.
 The tendency fields `Gⁿ` and `G⁻` can be specified via  optional `kwargs`.
 
-The scheme described by Le and Moin (1991) (see [LeMoin1991](@cite)). In a nutshel, the 3rd-order
+The scheme described by [LeMoin1991](@citet). In a nutshel, the 3rd-order
 Runge Kutta timestepper steps forward the state `Uⁿ` by `Δt` via 3 substeps. A pressure correction
 step is applied after at each substep.
 
diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl b/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl
index 8114d5b326..d49d7b7758 100644
--- a/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl
+++ b/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl
@@ -4,8 +4,8 @@ using Oceananigans.Operators
     AnisotropicMinimumDissipation{FT} <: AbstractTurbulenceClosure
 
 Parameters for the "anisotropic minimum dissipation" turbulence closure for large eddy simulation
-proposed originally by [Rozema15](@cite) and [Abkar16](@cite), then modified by [Verstappen18](@cite),
-and finally described and validated for by [Vreugdenhil18](@cite).
+proposed originally by [Rozema15](@citet) and [Abkar16](@citet), then modified by [Verstappen18](@citet),
+and finally described and validated for by [Vreugdenhil18](@citet).
 """
 struct AnisotropicMinimumDissipation{TD, PK, PN, PB} <: AbstractScalarDiffusivity{TD, ThreeDimensionalFormulation}
     Cν :: PN
@@ -58,9 +58,9 @@ Keyword arguments
         number or function is applied to all tracers. If a `NamedTuple`, it must possess
         a field specifying the Poncaré constant for every tracer.
 
-* `Cb`: Buoyancy modification multiplier (`Cb = nothing` turns it off, `Cb = 1` was used by [Abkar16](@cite)).
+* `Cb`: Buoyancy modification multiplier (`Cb = nothing` turns it off, `Cb = 1` was used by [Abkar16](@citet)).
         *Note*: that we _do not_ subtract the horizontally-average component before computing this
-        buoyancy modification term. This implementation differs from [Abkar16](@cite)'s proposal
+        buoyancy modification term. This implementation differs from [Abkar16](@citet)'s proposal
         and the impact of this approximation has not been tested or validated.
 
 By default: `C = Cν = Cκ = 1/12`, which is appropriate for a finite-volume method employing a