diff --git a/dapper/da_methods/baseline.py b/dapper/da_methods/baseline.py
index ed8196a5..87903539 100644
--- a/dapper/da_methods/baseline.py
+++ b/dapper/da_methods/baseline.py
@@ -1,6 +1,6 @@
 """Unsophisticated" but robust (widely applicable) DA methods.
 
-Many are based on `bib.raanes2016thesis`.
+Many are based on [raanes2016thesis][].
 """
 
 from typing import Callable, Optional
diff --git a/dapper/da_methods/ensemble.py b/dapper/da_methods/ensemble.py
index 1285e77e..96507316 100644
--- a/dapper/da_methods/ensemble.py
+++ b/dapper/da_methods/ensemble.py
@@ -28,7 +28,7 @@ class ens_method:
 class EnKF:
     """The ensemble Kalman filter.
 
-    Refs: `bib.evensen2009ensemble`.
+    Refs: [evensen2009ensemble][].
     """
 
     upd_a: str
@@ -67,8 +67,8 @@ def EnKF_analysis(E, Eo, hnoise, y, upd_a, stats=None, ko=None):
     This implementation includes several flavours and forms,
     specified by `upd_a`.
 
-    Main references: `bib.sakov2008deterministic`,
-    `bib.sakov2008implications`, `bib.hoteit2015mitigating`
+    Main references: [sakov2008deterministic][],
+    [sakov2008implications][], [hoteit2015mitigating][]
     """
     R = hnoise.C  # Obs noise cov
     N, Nx = E.shape  # Dimensionality
@@ -261,7 +261,7 @@ def post_process(E, infl, rot):
 def add_noise(E, dt, noise, method):
     """Treatment of additive noise for ensembles.
 
-    Refs: `bib.raanes2014ext`
+    Refs: [raanes2014ext][]
     """
     if noise.C == 0:
         return E
@@ -354,7 +354,7 @@ def sqrt_core():
 class EnKS:
     """The ensemble Kalman smoother.
 
-    Refs: `bib.evensen2009ensemble`
+    Refs: [evensen2009ensemble][]
 
     The only difference to the EnKF
     is the management of the lag and the reshapings.
@@ -414,7 +414,7 @@ def assimilate(self, HMM, xx, yy):
 class EnRTS:
     """EnRTS (Rauch-Tung-Striebel) smoother.
 
-    Refs: `bib.raanes2016thesis`
+    Refs: [raanes2016thesis][]
     """
 
     upd_a: str
@@ -484,7 +484,7 @@ def serial_inds(upd_a, y, cvR, A):
 class SL_EAKF:
     """Serial, covariance-localized EAKF.
 
-    Refs: `bib.karspeck2007experimental`.
+    Refs: [karspeck2007experimental][].
 
     In contrast with LETKF, this iterates over the observations rather
     than over the state (batches).
@@ -614,7 +614,7 @@ def local_analysis(ii):
 class LETKF:
     """Same as EnKF (Sqrt), but with localization.
 
-    Refs: `bib.hunt2007efficient`.
+    Refs: [hunt2007efficient][].
 
     NB: Multiproc. yields slow-down for [`mods.Lorenz96`][],
     even with `batch_size=(1,)`. But for [`mods.QG`][]
@@ -778,7 +778,7 @@ def hyperprior_coeffs(s, N, xN=1, g=0):
 
     - Reason 1: mode correction.
       These parameters bridge the Jeffreys (`xN=1`) and Dirac (`xN=Inf`) hyperpriors
-      for the prior covariance, B, as discussed in `bib.bocquet2015expanding`.
+      for the prior covariance, B, as discussed in [bocquet2015expanding][].
       Indeed, mode correction becomes necessary when $$ R \rightarrow \infty $$
       because then there should be no ensemble update (and also no inflation!).
       More specifically, the mode of `l1`'s should be adjusted towards 1
@@ -790,7 +790,7 @@ def hyperprior_coeffs(s, N, xN=1, g=0):
     - Reason 2: Boosting the inflation prior's certainty from N to xN*N.
       The aim is to take advantage of the fact that the ensemble may not
       have quite as much sampling error as a fully stochastic sample,
-      as illustrated in section 2.1 of `bib.raanes2019adaptive`.
+      as illustrated in section 2.1 of [raanes2019adaptive][].
 
     - Its damping effect is similar to work done by J. Anderson.
 
@@ -809,7 +809,7 @@ def hyperprior_coeffs(s, N, xN=1, g=0):
     eN = (N + 1) / N
     cL = (N + g) / N1
 
-    # Mode correction (almost) as in eqn 36 of `bib.bocquet2015expanding`
+    # Mode correction (almost) as in eqn 36 of [bocquet2015expanding][]
     prior_mode = eN / cL  # Mode of l1 (before correction)
     diagonal = pad0(s**2, N) + N1  # diag of Y@R.inv@Y + N1*I
     #                                           (Hessian of J)
@@ -847,7 +847,7 @@ def zeta_a(eN, cL, w):
 class EnKF_N:
     """Finite-size EnKF (EnKF-N).
 
-    Refs: `bib.bocquet2011ensemble`, `bib.bocquet2015expanding`
+    Refs: [bocquet2011ensemble][], [bocquet2015expanding][]
 
     This implementation is pedagogical, prioritizing the "dual" form.
     In consequence, the efficiency of the "primal" form suffers a bit.
@@ -985,7 +985,7 @@ def nvrs(w):
                 # l1 = 1.0
 
                 # Explicitly inflate prior
-                # => formulae look different from `bib.bocquet2015expanding`.
+                # => formulae look different from [bocquet2015expanding][].
                 A *= l1
                 Y *= l1
 
@@ -1000,7 +1000,7 @@ def nvrs(w):
                 else:
                     # Also include angular-radial co-dependence.
                     # Note: denominator not squared coz
-                    # unlike `bib.bocquet2015expanding` we have inflated Y.
+                    # unlike [bocquet2015expanding][] we have inflated Y.
                     Hw = (
                         Y @ R.inv @ Y.T / N1
                         + eye(N)
diff --git a/dapper/da_methods/other.py b/dapper/da_methods/other.py
index 95d44512..f13a997c 100644
--- a/dapper/da_methods/other.py
+++ b/dapper/da_methods/other.py
@@ -15,7 +15,7 @@
 class RHF:
     """Rank histogram filter.
 
-    Refs: `bib.anderson2010non`.
+    Refs: [anderson2010non][].
 
     Quick & dirty implementation without attention to (de)tails.
     """
@@ -86,9 +86,9 @@ def assimilate(self, HMM, xx, yy):
 class LNETF:
     """The Nonlinear-Ensemble-Transform-Filter (localized).
 
-    Refs: `bib.wiljes2016second`, `bib.todter2015second`.
+    Refs: [wiljes2016second][], [todter2015second][].
 
-    It is (supposedly) a deterministic upgrade of the NLEAF of `bib.lei2011moment`.
+    It is (supposedly) a deterministic upgrade of the NLEAF of [lei2011moment][].
     """
 
     N: int
diff --git a/dapper/da_methods/particle.py b/dapper/da_methods/particle.py
index 32ed799c..7e929329 100644
--- a/dapper/da_methods/particle.py
+++ b/dapper/da_methods/particle.py
@@ -23,7 +23,7 @@ class particle_method:
 class PartFilt:
     r"""Particle filter ≡ Sequential importance (re)sampling SIS (SIR).
 
-    Refs: `bib.wikle2007bayesian`, `bib.van2009particle`, `bib.chen2003bayesian`
+    Refs: [wikle2007bayesian][], [van2009particle][], [chen2003bayesian][]
 
     This is the bootstrap version: the proposal density is just
 
@@ -37,7 +37,7 @@ class PartFilt:
      - qroot: "Inflate" (anneal) the proposal noise kernels
        by this root to increase diversity.
        The weights are updated to maintain un-biased-ness.
-       See `bib.chen2003bayesian`, section VI-M.2
+       See [chen2003bayesian][], section VI-M.2
     """
 
     N: int
@@ -92,7 +92,7 @@ def assimilate(self, HMM, xx, yy):
 class OptPF:
     """'Optimal proposal' particle filter, also known as 'Implicit particle filter'.
 
-    Ref: `bib.bocquet2010beyond`.
+    Ref: [bocquet2010beyond][].
 
     !!! note Regularization (`Qs`) is here added BEFORE Bayes' rule.
         If `Qs==0`: OptPF should be equal to the bootstrap filter `PartFilt`.
@@ -520,7 +520,7 @@ def mask_unique_of_sorted(idx):
 def auto_bandw(N, M):
     """Optimal bandwidth (not bandwidth^2), as per Scott's rule-of-thumb.
 
-    Refs: `bib.doucet2001sequential` section 12.2.2, [Wik17]_ section "Rule_of_thumb"
+    Refs: [doucet2001sequential][] section 12.2.2, [Wik17]_ section "Rule_of_thumb"
     """
     return N ** (-1 / (M + 4))
 
@@ -534,7 +534,7 @@ def regularize(C12, E, idx, no_uniq_jitter):
     (so-named because Dirac-deltas are approximated  Gaussian kernels),
     which controls the strength of the jitter.
     This causes a bias. But, as N-->∞, the reg. bandwidth-->0, i.e. bias-->0.
-    Ref: `bib.doucet2001sequential`, section 12.2.2.
+    Ref: [doucet2001sequential][], section 12.2.2.
     """
     # Select
     E = E[idx]
@@ -554,7 +554,7 @@ def regularize(C12, E, idx, no_uniq_jitter):
 def resample(w, kind="Systematic", N=None, wroot=1.0):
     """Multinomial resampling.
 
-    Refs: `bib.doucet2009tutorial`, `bib.van2009particle`, `bib.liu2001theoretical`.
+    Refs: [doucet2009tutorial][], [van2009particle][], [liu2001theoretical][].
 
     - kind: 'Systematic', 'Residual' or 'Stochastic'.
       'Stochastic' corresponds to `rng.choice` or `rng.multinomial`.
@@ -571,7 +571,7 @@ def resample(w, kind="Systematic", N=None, wroot=1.0):
     - wroot: Adjust weights before resampling by this root to
       promote particle diversity and mitigate thinning.
       The outcomes of the resampling are then weighted to maintain un-biased-ness.
-      Ref: `bib.liu2001theoretical`, section 3.1
+      Ref: [liu2001theoretical][], section 3.1
 
     Note: (a) resampling methods are beneficial because they discard
     low-weight ("doomed") particles and reduce the variance of the weights.
diff --git a/dapper/da_methods/variational.py b/dapper/da_methods/variational.py
index 35b29045..9e344d17 100644
--- a/dapper/da_methods/variational.py
+++ b/dapper/da_methods/variational.py
@@ -28,10 +28,10 @@ class var_method:
 class iEnKS:
     """Iterative EnKS.
 
-    Special cases: EnRML, ES-MDA, iEnKF, EnKF `bib.raanes2019revising`.
+    Special cases: EnRML, ES-MDA, iEnKF, EnKF [raanes2019revising][].
 
-    As in `bib.bocquet2014iterative`, optimization uses Gauss-Newton.
-    See `bib.bocquet2012combining` for Levenberg-Marquardt.
+    As in [bocquet2014iterative][], optimization uses Gauss-Newton.
+    See [bocquet2012combining][] for Levenberg-Marquardt.
     If MDA=True, then there's not really any optimization,
     but rather Gaussian annealing.
 
@@ -41,13 +41,13 @@ class iEnKS:
 
         - "Sqrt"   : as in ETKF  , using a deterministic matrix square root transform.
         - "PertObs": as in EnRML , using stochastic, perturbed-observations.
-        - "Order1" : as in DEnKF of `bib.sakov2008deterministic`.
+        - "Order1" : as in DEnKF of [sakov2008deterministic][].
 
       Lag:
         Length of the DA window (DAW, multiples of dko, i.e. cycles).
 
-        - Lag=1 (default) => iterative "filter" iEnKF `bib.sakov2012iterative`.
-        - Lag=0           => maximum-likelihood filter `bib.zupanski2005maximum`.
+        - Lag=1 (default) => iterative "filter" iEnKF [sakov2012iterative][].
+        - Lag=0           => maximum-likelihood filter [zupanski2005maximum][].
 
       Shift : How far (in cycles) to slide the DAW.
               Fixed at 1 for code simplicity.
@@ -60,7 +60,7 @@ class iEnKS:
               Recommended: 1e-5.
 
       MDA   : Use iterations of the "multiple data assimlation" type.
-              Ref `bib.emerick2012history`
+              Ref [emerick2012history][]
 
       bundle: Use finite-diff. linearization instead of of least-squares regression.
               Makes the iEnKS very much alike the iterative, extended KF (IEKS).
@@ -70,8 +70,8 @@ class iEnKS:
     Total number of model simulations (of duration dto): N * (nIter*Lag + 1).
     (due to boundary cases: only asymptotically valid)
 
-    Refs: `bib.bocquet2012combining`, `bib.bocquet2013joint`,
-    `bib.bocquet2014iterative`.
+    Refs: [bocquet2012combining][], [bocquet2013joint][],
+    [bocquet2014iterative][].
     """
 
     upd_a: str
diff --git a/dapper/mods/KS/bocquet2019.py b/dapper/mods/KS/bocquet2019.py
index 8a285ca5..936ee162 100644
--- a/dapper/mods/KS/bocquet2019.py
+++ b/dapper/mods/KS/bocquet2019.py
@@ -1,4 +1,4 @@
-"""Settings as in `bib.bocquet2019consistency`."""
+"""Settings as in [bocquet2019consistency][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/LA/evensen2009.py b/dapper/mods/LA/evensen2009.py
index 651fa85c..26de4e63 100644
--- a/dapper/mods/LA/evensen2009.py
+++ b/dapper/mods/LA/evensen2009.py
@@ -1,4 +1,4 @@
-"""A mix of `bib.evensen2009ensemble` and `bib.sakov2008implications`.
+"""A mix of [evensen2009ensemble][] and [sakov2008implications][].
 
 !!! note
     Since there is no noise, and the system is stable,
diff --git a/dapper/mods/LA/raanes2015.py b/dapper/mods/LA/raanes2015.py
index b897cd03..b50835b5 100644
--- a/dapper/mods/LA/raanes2015.py
+++ b/dapper/mods/LA/raanes2015.py
@@ -1,4 +1,4 @@
-"""Reproduce results from Fig. 2 of `bib.raanes2014ext`."""
+"""Reproduce results from Fig. 2 of [raanes2014ext][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz05/__init__.py b/dapper/mods/Lorenz05/__init__.py
index 4d57e9f9..a6e9f15b 100644
--- a/dapper/mods/Lorenz05/__init__.py
+++ b/dapper/mods/Lorenz05/__init__.py
@@ -1,6 +1,6 @@
 """A multi-scale, smooth version of the classic Lorenz-96.
 
-This is an implementation of "Model III" of `bib.lorenz2005designing`.
+This is an implementation of "Model III" of [lorenz2005designing][].
 
 Similar to [`mods.LorenzUV`][] this model is designed
 to contain two different scales. However, in "Model III"
@@ -111,7 +111,7 @@ def step(self, E, t, dt):
 
 
 def summation_kernel(width):
-    """Prepare computation of the modified sum in `bib.lorenz2005designing`.
+    """Prepare computation of the modified sum in [lorenz2005designing][].
 
     Note: This gets repeatedly called, but actually the input is only ever
     `width = K` or `2*J`, so we should really cache or pre-compute.
@@ -131,7 +131,7 @@ def boxcar(x, n, method="direct"):
     For symmetry, if `n` is pair, the actual number of elements used is `n+1`,
     and the outer elements weighted by `0.5` to compensate for the `+1`.
 
-    This is the modified sum of `bib.lorenz2005designing`, used e.g. in eqn. 9.
+    This is the modified sum of [lorenz2005designing][], used e.g. in eqn. 9.
     For intuition: this type of summation (and the associated Sigma prime notation)
     may also be found for the "Trapezoidal rule" and in the inverse DFT used in
     spectral methods on a periodic domain.
@@ -216,7 +216,7 @@ def shift(x, k):
 
 
 def prodsum_self(x, k):
-    """Compute `prodsum(x, x, k)` efficiently: eqn 10 of `bib.lorenz2005designing`."""
+    """Compute `prodsum(x, x, k)` efficiently: eqn 10 of [lorenz2005designing][]."""
     W = boxcar(x, k)
     WW = shift(W, -2 * k) * shift(W, -k)
     WX = shift(W, -k) * shift(x, k)
diff --git a/dapper/mods/Lorenz63/anderson2010rhf.py b/dapper/mods/Lorenz63/anderson2010rhf.py
index dc05a8af..036cab44 100644
--- a/dapper/mods/Lorenz63/anderson2010rhf.py
+++ b/dapper/mods/Lorenz63/anderson2010rhf.py
@@ -1,4 +1,4 @@
-"""Settings from `bib.anderson2010non`."""
+"""Settings from [anderson2010non][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz63/bocquet2012.py b/dapper/mods/Lorenz63/bocquet2012.py
index 23c5f06e..a2641ee2 100644
--- a/dapper/mods/Lorenz63/bocquet2012.py
+++ b/dapper/mods/Lorenz63/bocquet2012.py
@@ -1,4 +1,4 @@
-"""Reproduce results from Fig 11 of `bib.bocquet2012combining`."""
+"""Reproduce results from Fig 11 of [bocquet2012combining][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz63/mandel2016.py b/dapper/mods/Lorenz63/mandel2016.py
index d9d74ddf..ce08a444 100644
--- a/dapper/mods/Lorenz63/mandel2016.py
+++ b/dapper/mods/Lorenz63/mandel2016.py
@@ -1,4 +1,4 @@
-"""Settings from `bib.mandel2016hybrid`."""
+"""Settings from [mandel2016hybrid][]."""
 
 import dapper.mods as modelling
 from dapper.mods.Lorenz63.sakov2012 import HMM as _HMM
diff --git a/dapper/mods/Lorenz63/sakov2012.py b/dapper/mods/Lorenz63/sakov2012.py
index 7dc1c699..11d8488d 100644
--- a/dapper/mods/Lorenz63/sakov2012.py
+++ b/dapper/mods/Lorenz63/sakov2012.py
@@ -1,4 +1,4 @@
-"""Reproduce results from Table 1 `bib.sakov2012iterative`."""
+"""Reproduce results from Table 1 [sakov2012iterative][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz63/wiljes2017.py b/dapper/mods/Lorenz63/wiljes2017.py
index d6def422..a1471de9 100644
--- a/dapper/mods/Lorenz63/wiljes2017.py
+++ b/dapper/mods/Lorenz63/wiljes2017.py
@@ -1,4 +1,4 @@
-"""Settings from `bib.wiljes2016second`."""
+"""Settings from [wiljes2016second][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz84/__init__.py b/dapper/mods/Lorenz84/__init__.py
index a069780f..4f27bf39 100644
--- a/dapper/mods/Lorenz84/__init__.py
+++ b/dapper/mods/Lorenz84/__init__.py
@@ -1,6 +1,6 @@
 """A chaotic system of size 3, like Lorenz-63, but with +complex geometry.
 
-Refs: `bib.lorenz1984irregularity`, `bib.lorenz2005look`
+Refs: [lorenz1984irregularity][], [lorenz2005look][]
 """
 
 import numpy as np
diff --git a/dapper/mods/Lorenz84/pajonk2012.py b/dapper/mods/Lorenz84/pajonk2012.py
index 35121dca..7b03b826 100644
--- a/dapper/mods/Lorenz84/pajonk2012.py
+++ b/dapper/mods/Lorenz84/pajonk2012.py
@@ -1,4 +1,4 @@
-"""Settings from `bib.pajonk2012deterministic`.
+"""Settings from [pajonk2012deterministic][].
 
 There is nothing to reproduce from the paper as there are no
 statistically converged numbers.
diff --git a/dapper/mods/Lorenz96/__init__.py b/dapper/mods/Lorenz96/__init__.py
index d0af3373..d80ef9d8 100644
--- a/dapper/mods/Lorenz96/__init__.py
+++ b/dapper/mods/Lorenz96/__init__.py
@@ -1,6 +1,6 @@
 """A 1D emulator of chaotic atmospheric behaviour.
 
-`bib.lorenz1996predictability`
+[lorenz1996predictability][]
 
 For a short introduction, see
 
diff --git a/dapper/mods/Lorenz96/anderson2009.py b/dapper/mods/Lorenz96/anderson2009.py
index e5bef2e7..8d8ea7e0 100644
--- a/dapper/mods/Lorenz96/anderson2009.py
+++ b/dapper/mods/Lorenz96/anderson2009.py
@@ -1,4 +1,4 @@
-"""A land-ocean setup from `bib.anderson2009spatially`."""
+"""A land-ocean setup from [anderson2009spatially][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz96/bocquet2010.py b/dapper/mods/Lorenz96/bocquet2010.py
index 70b5160d..d942e92b 100644
--- a/dapper/mods/Lorenz96/bocquet2010.py
+++ b/dapper/mods/Lorenz96/bocquet2010.py
@@ -1,4 +1,4 @@
-"""From Fig. 1 of `bib.bocquet2010beyond`."""
+"""From Fig. 1 of [bocquet2010beyond][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz96/bocquet2010_m40.py b/dapper/mods/Lorenz96/bocquet2010_m40.py
index 9727bb00..9909be9b 100644
--- a/dapper/mods/Lorenz96/bocquet2010_m40.py
+++ b/dapper/mods/Lorenz96/bocquet2010_m40.py
@@ -1,4 +1,4 @@
-"""From `bib.bocquet2010beyond` (again), but `ndim=40` (i.e. Fig. 5 of paper)."""
+"""From [bocquet2010beyond][] (again), but `ndim=40` (i.e. Fig. 5 of paper)."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz96/bocquet2015loc.py b/dapper/mods/Lorenz96/bocquet2015loc.py
index fed3abd6..866d773d 100644
--- a/dapper/mods/Lorenz96/bocquet2015loc.py
+++ b/dapper/mods/Lorenz96/bocquet2015loc.py
@@ -1,4 +1,4 @@
-"""Settings as in `bib.bocquet2016localization`."""
+"""Settings as in [bocquet2016localization][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz96/frei2013bridging.py b/dapper/mods/Lorenz96/frei2013bridging.py
index ddc5d1a4..f4c374f6 100644
--- a/dapper/mods/Lorenz96/frei2013bridging.py
+++ b/dapper/mods/Lorenz96/frei2013bridging.py
@@ -1,8 +1,8 @@
-"""Settings as in `bib.frei2013bridging`.
+"""Settings as in [frei2013bridging][].
 
 They also cite its use in the following:
 
-`bib.bengtsson2003toward`, `bib.lei2011moment`, `bib.frei2013mixture`.
+[bengtsson2003toward][], [lei2011moment][], [frei2013mixture][].
 """
 
 import numpy as np
diff --git a/dapper/mods/Lorenz96/hoteit2015.py b/dapper/mods/Lorenz96/hoteit2015.py
index 4fa0ac54..152b5084 100644
--- a/dapper/mods/Lorenz96/hoteit2015.py
+++ b/dapper/mods/Lorenz96/hoteit2015.py
@@ -1,4 +1,4 @@
-"""Reproduce results from `bib.hoteit2015mitigating`."""
+"""Reproduce results from [hoteit2015mitigating][]."""
 
 from dapper.mods.Lorenz96.sakov2008 import HMM as _HMM
 
diff --git a/dapper/mods/Lorenz96/miyoshi2011.py b/dapper/mods/Lorenz96/miyoshi2011.py
index 2c6763b6..ce41a7c6 100644
--- a/dapper/mods/Lorenz96/miyoshi2011.py
+++ b/dapper/mods/Lorenz96/miyoshi2011.py
@@ -1,6 +1,6 @@
 """A land-ocean setup.
 
-Refs: `bib.miyoshi2011gaussian`, which was inspired by `bib.lorenz1998optimal`.
+Refs: [miyoshi2011gaussian][], which was inspired by [lorenz1998optimal][].
 """
 
 import numpy as np
diff --git a/dapper/mods/Lorenz96/pinheiro2019.py b/dapper/mods/Lorenz96/pinheiro2019.py
index 55074c4b..c8c924d0 100644
--- a/dapper/mods/Lorenz96/pinheiro2019.py
+++ b/dapper/mods/Lorenz96/pinheiro2019.py
@@ -1,4 +1,4 @@
-"""Settings from `bib.pinheiro2019efficient`."""
+"""Settings from [pinheiro2019efficient][]."""
 
 import dapper.mods as modelling
 import dapper.mods.Lorenz96 as model
diff --git a/dapper/mods/Lorenz96/raanes2016.py b/dapper/mods/Lorenz96/raanes2016.py
index 427bcb06..95b185db 100644
--- a/dapper/mods/Lorenz96/raanes2016.py
+++ b/dapper/mods/Lorenz96/raanes2016.py
@@ -1,4 +1,4 @@
-"""Reproduce `bib.raanes2015rts`."""
+"""Reproduce [raanes2015rts][]."""
 
 import dapper.mods as modelling
 from dapper.mods.Lorenz96.sakov2008 import HMM as _HMM
diff --git a/dapper/mods/Lorenz96/sakov2008.py b/dapper/mods/Lorenz96/sakov2008.py
index f28f6373..58f3c795 100644
--- a/dapper/mods/Lorenz96/sakov2008.py
+++ b/dapper/mods/Lorenz96/sakov2008.py
@@ -1,9 +1,9 @@
-"""Settings as in `bib.sakov2008deterministic`.
+"""Settings as in [sakov2008deterministic][].
 
 This HMM is used (with small variations) in many DA papers, for example
 
-`bib.bocquet2011ensemble`, `bib.sakov2012iterative`,
-`bib.bocquet2015expanding`, `bib.bocquet2013joint`.
+[bocquet2011ensemble][], [sakov2012iterative][],
+[bocquet2015expanding][], [bocquet2013joint][].
 """
 
 import numpy as np
diff --git a/dapper/mods/Lorenz96/todter2015.py b/dapper/mods/Lorenz96/todter2015.py
index e3bc35f5..522069fa 100644
--- a/dapper/mods/Lorenz96/todter2015.py
+++ b/dapper/mods/Lorenz96/todter2015.py
@@ -1,4 +1,4 @@
-"""Concerns figure 4 of `bib.todter2015second`."""
+"""Concerns figure 4 of [todter2015second][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/Lorenz96s/__init__.py b/dapper/mods/Lorenz96s/__init__.py
index b99a620d..5351b45c 100644
--- a/dapper/mods/Lorenz96s/__init__.py
+++ b/dapper/mods/Lorenz96s/__init__.py
@@ -1,6 +1,6 @@
 """A perfect-random version of Lorenz-96.
 
-Used by `bib.grudzien2020numerical` to study the precision of
+Used by [grudzien2020numerical][] to study the precision of
 stochastic integration schemes.
 
 Both the model and truth are to be integrated by the same *random* model (with almost
@@ -11,7 +11,7 @@
 
 The truth twin should be generated by the order 2.0 Taylor scheme below,
 for the accuracy with respect to convergence in the strong sense.
-See `bib.grudzien2020numerical` for a full discussion of benchmarks on this
+See [grudzien2020numerical][] for a full discussion of benchmarks on this
 model and statistically robust configurations.
 
 This study uses no multiplicative inflation / localization or other
@@ -19,7 +19,7 @@
 observation EnKF as a simple estimator to study the asymptotic filtering
 statistics under different model scenarios.
 
-The purpose of the study in `bib.grudzien2020numerical` was to explore the relationships
+The purpose of the study in [grudzien2020numerical][] was to explore the relationships
 between:
 
 - numerical discretization error in truth twins;
@@ -38,7 +38,7 @@
 thus corresponds to a wider variance of the relizations of the diffeomorphsims
 that generate the model / truth twin between observation times.
 
-It is demonstrated by `bib.grudzien2020numerical` that the model error due to
+It is demonstrated by [grudzien2020numerical][] that the model error due to
 discretization of the SDE equations of motion is most detrimental to the filtering cycle
 when model uncertainty is low and observation precision is high.  In other
 configurations, such as those with high model uncertainty, the differences between
@@ -93,7 +93,7 @@ def l96s_tay2_step(x, t, dt, s):
     model twin will be generated by on of the wrappers below.  The order 2.0
     Taylor-Stratonovich discretization scheme implemented here is the basic formulation
     which makes a Fourier truncation at p=1 for the Brownian bridge process. See
-    `bib.grudzien2020numerical` for full details of the scheme and other versions."""
+    [grudzien2020numerical][] for full details of the scheme and other versions."""
 
     # Infer system dimension
     sys_dim = len(x)
diff --git a/dapper/mods/Lorenz96s/grudzien2020.py b/dapper/mods/Lorenz96s/grudzien2020.py
index 50da5464..d6dcd2bc 100644
--- a/dapper/mods/Lorenz96s/grudzien2020.py
+++ b/dapper/mods/Lorenz96s/grudzien2020.py
@@ -1,4 +1,4 @@
-"""Settings as in `bib.grudzien2020numerical`."""
+"""Settings as in [grudzien2020numerical][]."""
 
 import dapper.mods as modelling
 from dapper.mods.Lorenz96 import Tplot, x0
diff --git a/dapper/mods/LorenzUV/__init__.py b/dapper/mods/LorenzUV/__init__.py
index c2aba4f6..30f54245 100644
--- a/dapper/mods/LorenzUV/__init__.py
+++ b/dapper/mods/LorenzUV/__init__.py
@@ -1,6 +1,6 @@
 """The 2-scale/layer/speed coupled version of Lorenz-96.
 
-See `bib.wilks2005effects`
+See [wilks2005effects][]
 - U:  large amp, low frequency vars: convective events
 - V:  small amp, high frequency vars: large-scale synoptic events
 
@@ -45,7 +45,7 @@ def newfun(x, *args, reverse=False, **kwargs):
 class model_instance:
     """Use OOP to facilitate having multiple parameter settings simultaneously.
 
-    Default parameters from `bib.wilks2005effects`.
+    Default parameters from [wilks2005effects][].
     """
 
     def __init__(self, nU=8, J=32, F=20, h=1, b=10, c=10):
diff --git a/dapper/mods/LorenzUV/lorenz96.py b/dapper/mods/LorenzUV/lorenz96.py
index aa119ea2..6b75a315 100644
--- a/dapper/mods/LorenzUV/lorenz96.py
+++ b/dapper/mods/LorenzUV/lorenz96.py
@@ -1,4 +1,4 @@
-"""As in `bib.lorenz1996predictability`."""
+"""As in [lorenz1996predictability][]."""
 
 from pathlib import Path
 
diff --git a/dapper/mods/LorenzUV/wilks05.py b/dapper/mods/LorenzUV/wilks05.py
index fa347bbd..43a671c3 100644
--- a/dapper/mods/LorenzUV/wilks05.py
+++ b/dapper/mods/LorenzUV/wilks05.py
@@ -1,4 +1,4 @@
-"""Uses `nU`, `J`, `F` as in [`mods.LorenzUV`][] ie. from `bib.wilks2005effects`.
+"""Uses `nU`, `J`, `F` as in [`mods.LorenzUV`][] ie. from [wilks2005effects][].
 
 Obs settings taken from different places (=> quasi-linear regime).
 """
diff --git a/dapper/mods/LotkaVolterra/__init__.py b/dapper/mods/LotkaVolterra/__init__.py
index fe83f1d5..42f61a82 100644
--- a/dapper/mods/LotkaVolterra/__init__.py
+++ b/dapper/mods/LotkaVolterra/__init__.py
@@ -2,7 +2,7 @@
 
 Refs:
 [Wiki](https://en.wikipedia.org/wiki/Competitive_Lotka-Volterra_equations),
-`bib.vano2006chaos`.
+[vano2006chaos][].
 """
 
 import numpy as np
diff --git a/dapper/mods/QG/__init__.py b/dapper/mods/QG/__init__.py
index ee6dc5e7..50c73c09 100644
--- a/dapper/mods/QG/__init__.py
+++ b/dapper/mods/QG/__init__.py
@@ -1,4 +1,4 @@
-"""Quasi-geostraphic 2D flow. Described in detail by `bib.sakov2008deterministic`.
+"""Quasi-geostraphic 2D flow. Described in detail by [sakov2008deterministic][].
 
 Adapted from Pavel Sakov's enkf-matlab package.
 
diff --git a/dapper/mods/QG/counillon2009.py b/dapper/mods/QG/counillon2009.py
index 1d6f6436..aac81145 100644
--- a/dapper/mods/QG/counillon2009.py
+++ b/dapper/mods/QG/counillon2009.py
@@ -1,4 +1,4 @@
-"""Reproduce experiments from `bib.counillon2009application`."""
+"""Reproduce experiments from [counillon2009application][]."""
 
 import dapper.mods as modelling
 from dapper.mods.QG import model_config
diff --git a/dapper/mods/QG/sakov2008.py b/dapper/mods/QG/sakov2008.py
index 604575f1..fc45f5e0 100644
--- a/dapper/mods/QG/sakov2008.py
+++ b/dapper/mods/QG/sakov2008.py
@@ -1,4 +1,4 @@
-"""Reproduce results from `bib.sakov2008deterministic`."""
+"""Reproduce results from [sakov2008deterministic][]."""
 
 import numpy as np
 
diff --git a/dapper/mods/VL20/__init__.py b/dapper/mods/VL20/__init__.py
index 138093ab..cdbe2544 100644
--- a/dapper/mods/VL20/__init__.py
+++ b/dapper/mods/VL20/__init__.py
@@ -1,7 +1,7 @@
 """
 Single-scale Lorenz-96 with an added thermodynamic component.
 
-Refs: `bib.vissio2020mechanics`
+Refs: [vissio2020mechanics][]
 """
 
 import numpy as np
diff --git a/dapper/mods/VL20/demo.py b/dapper/mods/VL20/demo.py
index 9885d30c..55ec3d28 100644
--- a/dapper/mods/VL20/demo.py
+++ b/dapper/mods/VL20/demo.py
@@ -1,6 +1,6 @@
 """Demonstrate the Vissio-Lucarini-20 model.
 
-Reproduce Hovmoller diagram Fig 4. in `bib.vissio2020mechanics`.
+Reproduce Hovmoller diagram Fig 4. in [vissio2020mechanics][].
 """
 
 import numpy as np
diff --git a/dapper/mods/integration.py b/dapper/mods/integration.py
index 8c82c2ec..f34f16de 100644
--- a/dapper/mods/integration.py
+++ b/dapper/mods/integration.py
@@ -19,7 +19,7 @@ def rk4(f, x, t, dt, stages=4, s=0):
     For SDEs with additive noise (`s>0`), the order of convergence
     (both weak and strong) is 1 for `stages` equal to 1 or 4.
     These correspond to the classic Euler-Maruyama scheme and the Runge-Kutta
-    scheme for S-ODEs respectively, see `bib.grudzien2020numerical`
+    scheme for S-ODEs respectively, see [grudzien2020numerical][]
     for a DA-specific discussion on integration schemes and their discretization errors.
 
     Parameters
diff --git a/docs/examples/param_estim.py b/docs/examples/param_estim.py
index 85184e73..99dc59c5 100644
--- a/docs/examples/param_estim.py
+++ b/docs/examples/param_estim.py
@@ -6,7 +6,7 @@
 # (concatenating) the state vector with the parameter (vector).
 #
 # This particular experiment is a reproduction of section 3.3 of
-# `bib.bocquet2013joint`. Briefly: the unknown parameter is the forcing used in
+# [bocquet2013joint][]. Briefly: the unknown parameter is the forcing used in
 # Lorenz-96; only the state is observed; both parameter and state get estimated.
 #
 # This example builds mostly on `examples/basic_2.py`. For brevity, it does not
diff --git a/docs/examples/stoch_model1.py b/docs/examples/stoch_model1.py
index 5ff1c7ce..95156a78 100644
--- a/docs/examples/stoch_model1.py
+++ b/docs/examples/stoch_model1.py
@@ -1,4 +1,4 @@
-# ## Reproduce basics of `bib.grudzien2020numerical`, and do live plotting.
+# ## Reproduce basics of [grudzien2020numerical][], and do live plotting.
 # This is a simple demonstration of the effect of using low-precision numerics
 # for a perfect-random model configuration in which the system is not dominated
 # by the noise.  In this case, we see a failure of the ensemble generated