diff --git a/dapper/mods/NS2D/2D_Periodic_NavierStokes_Solver_for_Data_Assimilation_with_DAPPER.pdf b/dapper/mods/NS2D/2D_Periodic_NavierStokes_Solver_for_Data_Assimilation_with_DAPPER.pdf
new file mode 100644
index 00000000..7d5eca7b
Binary files /dev/null and b/dapper/mods/NS2D/2D_Periodic_NavierStokes_Solver_for_Data_Assimilation_with_DAPPER.pdf differ
diff --git a/dapper/mods/NS2D/__init__.py b/dapper/mods/NS2D/__init__.py
new file mode 100644
index 00000000..15a254cb
--- /dev/null
+++ b/dapper/mods/NS2D/__init__.py
@@ -0,0 +1,317 @@
+"""The Navier-Stokes equations in 2D. Using the Streamfunction form of these equations:
+laplacian(psi)_t = psi_y * laplacian(psi)_x - psi_x * laplacian(psi)_y +
+nu * laplacian(laplacian(psi))
+
+See demo.py for an example of a Taylor-Green vortex"""
+
+import matplotlib as mpl
+import numpy as np
+from scipy.differentiate import jacobian
+
+import dapper.tools.liveplotting as LP
+from dapper.dpr_config import DotDict
+
+N_global = 64  # This is awkward but necessary because otherwise spatial2D would have
+# to be modified. If there is a better way to do this please change
+
+
+def Model(N=64, Lxy=2 * np.pi, dt=0.0001, nu=0.01, T=0.1):
+    """
+    Parameters
+    ----------
+    N : int
+        N = Nx* = Ny* = number of grid points in each dimension.
+    Lxy : float
+        Domain size in each dimension.
+    dt : float
+        Time step.
+    nu : float
+        Kinematic viscosity (dimensionless units); inv. prop. to Reynolds number.
+    T : float
+        Experiment time; total sim steps = T/dt.
+
+    Returns
+    -------
+    dd : DotDict
+        Model parameters and functions.
+
+    Notes
+    -----
+    See the report pdf in the folder for further explanation of the implementation
+    and mathematics\n
+    *Nx and Ny are not the same as the convention in DAPPER where Nx is the state size
+    and Ny is the observation size; here they just represent N. This is also the
+    convention in the report. More confusingly, sometimes N represents the grid size
+    and sometimes represents ensemble size; each function will delineate what N is
+    in its parameters if present.
+    """
+    assert (
+        N == N_global
+    ), f"Warning: change N_global in __init__.py to match your parameter ({N})"
+
+    h = dt  # alias -- prevents o/w in step()
+
+    # Initialize grid and wavenumbers (N = number of grid points in each dimension)
+    x = np.linspace(0, Lxy, N, False)
+    y = np.linspace(0, Lxy, N, False)
+    X, Y = np.meshgrid(x, y)
+
+    # Wavenumbers
+    kk = np.fft.fftfreq(N, Lxy / N) * 2 * np.pi
+    kx = kk.reshape((N, 1))
+    ky = kk.reshape((1, N))
+    k2 = kx**2 + ky**2
+    k2[0, 0] = 1  # avoid division by zero
+    # initial conditions for a Taylor-Green vortex; to use decaying Kolmogorov,
+    # comment 56 and uncomment 57-59
+    psi = np.sin(X) * np.sin(Y)
+    # U = 1
+    # k = 2
+    # psi = -U/k * np.cos(k*Y)
+    x0 = psi.copy()
+    L = -nu * k2  # generalize
+    L_flat = L.flatten()
+    E = np.exp(h * L_flat)
+    E2 = np.exp(h * L_flat / 2)
+
+    E = E.reshape((N, N))
+    E2 = E2.reshape((N, N))
+
+    def dealias(f_hat):
+        """Implements the 2/3 rule for dealiasing (Orszag, 1971,
+        https://doi.org/10.1175/1520-0469(1971)028%3C1074:OTEOAI%3E2.0.CO;2)\n
+        See section 2.4 of report
+
+        Parameters
+        ----------
+        f_hat : ndarray
+            2D array in frequency space.
+
+        Returns
+        -------
+        ndarray
+            Dealiased 2D array in frequency space.
+
+        """
+        N = f_hat.shape[0]
+        k_cutoff = N // 3
+
+        k_int = np.fft.fftfreq(N) * N
+        mask_1d = np.abs(k_int) < k_cutoff
+        mask_2d = np.outer(mask_1d, mask_1d)
+        return f_hat * mask_2d  # element-wise multiplication
+
+    # J = dpsi/dy * domega/dx - dpsi/dx * domega/dy
+    def det_jacobian_equation(psi_hat):
+        """
+        Computes J in the streamfunction formulation of the 2D NSE \n
+        See section 2.1, equation 5 of report
+
+        Parameters
+        ----------
+        psi_hat : ndarray
+            2D array in frequency space.
+
+        Returns
+        -------
+        ndarray
+            J in physical space.
+        """
+        scale = 1 / (Lxy * Lxy)  # Rescale to match the original domain size
+        dpsi_dy = np.fft.ifft2(1j * ky * psi_hat) * scale
+        domega_dx = np.fft.ifft2(1j * kx * k2 * psi_hat) * scale
+        dpsi_dx = np.fft.ifft2(1j * kx * psi_hat) * scale
+        domega_dy = np.fft.ifft2(1j * ky * k2 * psi_hat) * scale
+
+        return dpsi_dy * domega_dx - dpsi_dx * domega_dy
+
+    # ETD-RK4 method used in the KS model; see dapper/mods/KS/__init__.py
+    # Based on kursiv.m of Kassam and Trefethen, 2002,
+    # doi.org/10.1137/S1064827502410633.
+
+    # Adapted for the Navier-Stokes equations in 2D.
+    def NL(psi_hat):
+        """Returns the nonlinear term `N` in equation 11 of report
+        Parameters
+        ----------
+        psi_hat : ndarray
+            2D array in frequency space.
+
+        Returns
+        -------
+        ndarray
+            Nonlinear term in frequency space.
+        """
+        J = det_jacobian_equation(dealias(psi_hat))
+        J_hat = np.fft.fft2(J)
+        return -J_hat
+
+    # For the step function
+    def f(psi):  # consider redefining with psi_hat
+        return np.fft.ifft2(NL(np.fft.fft2(psi)) + nu * k2 * k2 * np.fft.fft2(psi)).real
+
+    def dstep_dx(psi, t, dt):
+        return jacobian(f, psi)
+
+    # Contour integral approximation and coefficients
+    nRoots = 16
+    roots = np.exp(1j * np.pi * (0.5 + np.arange(1, nRoots + 1)) / nRoots)
+
+    CL = h * L_flat[:, None] + roots
+    Q = h * ((np.exp(CL / 2) - 1) / CL).mean(axis=-1).real
+
+    f1 = h * ((-4 - CL + np.exp(CL) * (4 - 3 * CL + CL**2)) / CL**3).mean(axis=-1).real
+    f2 = h * ((2 + CL + np.exp(CL) * (-2 + CL)) / CL**3).mean(axis=-1).real
+    f3 = h * ((-4 - 3 * CL - CL**2 + np.exp(CL) * (4 - CL)) / CL**3).mean(axis=-1).real
+
+    Q = Q.reshape((1, N, N))
+    f1 = f1.reshape((1, N, N))
+    f2 = f2.reshape((1, N, N))
+    f3 = f3.reshape((1, N, N))
+
+    def step_ETD_RK4(x, t, dt):
+        """
+        Step function using ETD-RK4; x = psi
+
+        Parameters
+        ----------
+        x : ndarray
+            flattened 2D array in physical space.
+        t : float
+        dt : float
+            time step; must match initialized dt.
+
+        Returns
+        -------
+        ndarray
+            flattened 2D array in physical space after time step."""
+        epsilon = 1e-6
+        assert abs(dt - h) < epsilon, "dt must match the initialized dt"
+
+        psi = x.reshape((N, N))
+        p_hat = np.fft.fft2(psi)
+        # omega = laplacian(psi)
+        w_hat = k2 * p_hat
+        N1 = NL(p_hat)
+        v1 = E2 * w_hat + Q * N1
+        N2a = NL(v1)
+        v2a = E2 * w_hat + Q * N2a
+        N2b = NL(v2a)
+        v2b = E2 * v1 + Q * (2 * N2b - N1)
+        N3 = NL(v2b)
+        omega_hat_new = E * w_hat + N1 * f1 + 2 * (N2a + N2b) * f2 + N3 * f3
+        psi_hat_new = omega_hat_new / k2  # i^2 = -1; -1 / - 1 = 1.
+        psi_hat_new[0, 0] = 0  # Enforce zero mean for
+        psi_new = np.fft.ifft2(psi_hat_new).real
+        return psi_new.flatten()
+
+    # Vectorized versions of above functions for ensemble computations
+    def det_jacobian_equation_vec(psi_hat_batch):
+        # psi_hat_batch: (N_ens, N, N)
+        scale = 1 / (Lxy * Lxy)
+        dpsi_dy = np.fft.ifft2(1j * ky * psi_hat_batch, axes=(-2, -1)) * scale
+        domega_dx = np.fft.ifft2(1j * kx * k2 * psi_hat_batch, axes=(-2, -1)) * scale
+        dpsi_dx = np.fft.ifft2(1j * kx * psi_hat_batch, axes=(-2, -1)) * scale
+        domega_dy = np.fft.ifft2(1j * ky * k2 * psi_hat_batch, axes=(-2, -1)) * scale
+        return dpsi_dy * domega_dx - dpsi_dx * domega_dy
+
+    def dealias_vec(f_hat_batch):
+        # f_hat_batch: (N_ens, N, N)
+        Nf = f_hat_batch.shape[-1]
+        k_cutoff = Nf // 3
+        k_int = np.fft.fftfreq(Nf) * Nf
+        mask_1d = np.abs(k_int) < k_cutoff
+        mask_2d = np.outer(mask_1d, mask_1d)
+        return f_hat_batch * mask_2d  # broadcasting
+
+    def nonlinear_vec(psi_hat):
+        # psi_hat: (N_ens, N, N)
+        J = det_jacobian_equation_vec(dealias_vec(psi_hat))
+        J_hat = np.fft.fft2(J, axes=(-2, -1))
+        return -J_hat
+
+    def step_ETD_RK4_vec(X, t, dt):
+        # X: (N_ens, N*N) or (N_ens, N, N)
+        if X.ndim == 2 and X.shape[1] == N * N:
+            X = X.reshape((-1, N, N))
+        N_ens = X.shape[0]
+        p_hat = np.fft.fft2(X, axes=(-2, -1))
+        # omega = laplacian(psi)
+        w_hat = k2 * p_hat
+        N1 = nonlinear_vec(p_hat)
+        v1 = E2 * w_hat + Q * N1
+        N2a = nonlinear_vec(v1)
+        v2a = E2 * w_hat + Q * N2a
+        N2b = nonlinear_vec(v2a)
+        v2b = E2 * v1 + Q * (2 * N2b - N1)
+        N3 = nonlinear_vec(v2b)
+        omega_hat_new = E * w_hat + N1 * f1 + 2 * (N2a + N2b) * f2 + N3 * f3
+        psi_hat_new = omega_hat_new / k2  # i^2 = -1; -1 / - 1 = 1.
+        psi_hat_new[0, 0] = 0  # Enforce zero mean for
+        psi_new = np.fft.ifft2(psi_hat_new).real
+        return psi_new.reshape((N_ens, N * N))
+
+    def step_parallel(E, t, dt):
+        """Parallelized step for ensemble (2D array); otherwise normal step for
+        single state (1D).
+        Parameters
+        ----------
+        E : ndarray
+            Either 1D single state or 2D ensemble of states.
+        t : float
+        dt : float
+
+        Returns
+        -------
+        ndarray
+            Ensemble or single state advanced one time step."""
+        if E.ndim == 1:
+            return step_ETD_RK4(E, t, dt)
+        if E.ndim == 2:
+            return step_ETD_RK4_vec(E, t, dt)
+
+    dd = DotDict(
+        dt=dt,
+        nu=nu,
+        DL=2,
+        step=step_parallel,  # Use the parallelized step function when possible
+        dstep_dx=dstep_dx,
+        Nx=N,
+        x0=x0,
+        T=T,
+    )
+    return dd
+
+
+##Liveplotting mostly copied from QG model
+
+
+def square(x):
+    return x.reshape(N_global, N_global)
+
+
+def ind2sub(ind):
+    return np.unravel_index(ind, (N_global, N_global))
+
+
+cm = mpl.cm.viridis
+center = N_global * int(N_global / 2) + int(0.5 * N_global)
+
+
+def LP_setup(jj):
+    return [
+        (
+            1,
+            LP.spatial2d(
+                square,
+                ind2sub,
+                jj,
+                cm,
+                clims=((-1, 1), (-1, 1), (-1, 1), (-1, 1)),
+                domainlims=(2 * np.pi, 2 * np.pi),
+                periodic=(True, True),
+            ),
+        ),
+        (0, LP.spectral_errors),
+        (0, LP.sliding_marginals),
+    ]
diff --git a/dapper/mods/NS2D/demo.py b/dapper/mods/NS2D/demo.py
new file mode 100644
index 00000000..35f4baab
--- /dev/null
+++ b/dapper/mods/NS2D/demo.py
@@ -0,0 +1,47 @@
+import matplotlib.animation as animation
+import matplotlib.pyplot as plt
+import numpy as np
+import tqdm
+
+from dapper.mods.NS2D import Model
+
+model = Model(N=64, dt=0.01, T=200, nu=1 / 1600)
+T = model.T
+tt = np.linspace(0, T, int(T / model.dt), endpoint=True, dtype=np.float64)
+EE = np.zeros((len(tt), model.Nx * model.Nx))
+EE[0] = model.x0.flatten()  # IC comes from model; change IC to change demo output
+
+for k in tqdm.tqdm(range(1, len(tt))):
+    EE[k] = model.step(EE[k - 1], np.nan, model.dt)
+
+
+def animate_snapshots(psis, snapshot_steps):
+    # Select n evenly spaced indices (n = snapshot_steps)
+    indices = np.linspace(0, len(psis), snapshot_steps, False, dtype=int)
+    psi_snapshots = psis[indices]
+    # For title, get the actual time step for each snapshot
+    step_numbers = indices
+
+    # Animate
+    fig, ax = plt.subplots()
+    im = ax.imshow(
+        psi_snapshots[0],
+        origin="lower",
+        cmap="viridis",
+        extent=(0, model.DL * np.pi, 0, model.DL * np.pi),
+    )
+    ax.set_title(f"Streamfunction ψ, step {step_numbers[0]}")
+    plt.colorbar(im, ax=ax)
+
+    def update(frame):
+        im.set_data(psi_snapshots[frame])
+        ax.set_title(f"Streamfunction ψ, step {step_numbers[frame]}")
+        return (im,)
+
+    _ = animation.FuncAnimation(
+        fig, update, frames=len(psi_snapshots), interval=100, blit=False
+    )
+    plt.show()
+
+
+animate_snapshots(EE.reshape(len(tt), model.Nx, model.Nx), 10)
diff --git a/dapper/mods/NS2D/example_paper.py b/dapper/mods/NS2D/example_paper.py
new file mode 100644
index 00000000..a9a0c09d
--- /dev/null
+++ b/dapper/mods/NS2D/example_paper.py
@@ -0,0 +1,85 @@
+# Just using bocquet2019 settings for format and example;
+# reproduce paper results if found
+
+import numpy as np
+
+import dapper.mods as modelling
+from dapper.mods.NS2D import LP_setup, Model
+from dapper.tools.localization import nd_Id_localization
+
+noise_amp = 0.005
+System = Model(T=0.1, N=64, dt=0.0001, nu=0.01)
+Nx = System.Nx
+
+tseq = modelling.Chronology(System.dt, dko=1, BurnIn=0.01, T=0.1)
+
+Dyn = {
+    "M": np.prod((Nx, Nx)),
+    "model": System.step,
+    "linear": System.dstep_dx,
+    "noise": 0,
+}
+X0 = modelling.RV(
+    M=Dyn["M"],
+    func=lambda N: System.x0.flatten()[None, :]
+    + noise_amp
+    * np.random.randn(
+        N, Dyn["M"]
+    ),  # Initial perturbation is noise_amp * N(0, 1) for each gridpoint
+)
+
+
+Obs = modelling.Id_Obs(Nx**2)
+Obs["noise"] = 1
+
+Obs["localizer"] = nd_Id_localization((Nx,), (4,))
+
+rstream = np.random.RandomState()
+jj = modelling.linspace_int(Nx, Nx)
+max_offset = jj[1] - jj[0]
+
+
+def obs_inds(ko):
+    def random_offset():
+        rstream.seed(ko)
+        u = rstream.rand()
+        return int(np.floor(max_offset * u))
+
+    return jj + random_offset()
+
+
+def obs_now(ko):
+    jj = obs_inds(ko)
+    shape = (Nx, Nx)
+
+    @modelling.ens_compatible
+    def hmod(E):
+        return E[jj]
+
+    # Localization.
+    batch_shape = [4, 4]  # width (in grid points) of each state batch.
+    # Increasing the width
+    #  => quicker analysis (but less rel. speed-up by parallelzt., depending on NPROC)
+    #  => worse (increased) rmse (but width 4 is only slightly worse than 1);
+    #     if inflation is applied locally, then rmse might actually improve.
+    localizer = nd_Id_localization((shape)[::-1], batch_shape[::-1], jj, periodic=False)
+
+    Obs = {
+        "M": Nx,
+        "model": hmod,
+        "noise": modelling.GaussRV(C=4 * np.eye(Nx)),
+        "localizer": localizer,
+    }
+
+    # Moving localization mask for smoothers:
+    Obs["loc_shift"] = lambda ii, dt: ii  # no movement (suboptimal, but easy)
+
+    # Jacobian left unspecified coz it's (usually) employed by methods that
+    # compute full cov, which in this case is too big.
+
+    return modelling.Operator(**Obs)
+
+
+Obs = dict(time_dependent=lambda ko: obs_now(ko))
+HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0, LP=LP_setup(obs_inds))
+# HMM.liveplotters = LP_setup(obs_inds)
diff --git a/dapper/tools/liveplotting.py b/dapper/tools/liveplotting.py
index 753ad53f..7f9bad1e 100644
--- a/dapper/tools/liveplotting.py
+++ b/dapper/tools/liveplotting.py
@@ -1262,6 +1262,8 @@ def spatial2d(
     obs_inds=(),
     cm=plt.cm.jet,
     clims=((-40, 40), (-40, 40), (-10, 10), (-10, 10)),
+    domainlims=(2 * np.pi, 2 * np.pi),  # Lx, Ly
+    periodic=(True, True),
 ):
     def init(fignum, stats, key0, plot_u, E, P, **kwargs):
         GS = {"left": 0.125 - 0.04, "right": 0.9 - 0.04}
@@ -1301,11 +1303,43 @@ def init(fignum, stats, key0, plot_u, E, P, **kwargs):
 
         # Plot
         # - origin='lower' might get overturned by set_ylim() below.
-        im_11 = ax_11.imshow(square(mu[key0]), cmap=cm)
-        im_12 = ax_12.imshow(square(xx[k]), cmap=cm)
-        # hot is better, but needs +1 colorbar
-        im_21 = ax_21.imshow(square(spread[key0]), cmap=plt.cm.bwr)
-        im_22 = ax_22.imshow(square(err[key0]), cmap=plt.cm.bwr)
+        arr_mu = square(mu[key0])
+        arr_xx = square(xx[k])
+        im_11 = ax_11.imshow(
+            arr_mu,
+            cmap=cm,
+            extent=(0, domainlims[0], 0, domainlims[1]),
+            origin="lower",
+            interpolation="nearest",
+        )
+        im_12 = ax_12.imshow(
+            arr_xx,
+            cmap=cm,
+            extent=(0, domainlims[0], 0, domainlims[1]),
+            origin="lower",
+            interpolation="nearest",
+        )
+        # Ensure spread and error arrays are reshaped correctly
+        arr_spread = square(spread[key0])
+        arr_err = square(err[key0])
+
+        # Update imshow for spread and error plots
+        im_21 = ax_21.imshow(
+            arr_spread,
+            cmap=plt.cm.bwr,
+            extent=(0, domainlims[0], 0, domainlims[1]),
+            origin="lower",
+            interpolation="nearest",
+        )
+        im_22 = ax_22.imshow(
+            arr_err,
+            cmap=plt.cm.bwr,
+            extent=(0, domainlims[0], 0, domainlims[1]),
+            origin="lower",
+            interpolation="nearest",
+        )
+
+        # Update ims tuple
         ims = (im_11, im_12, im_21, im_22)
         # Obs init -- a list where item 0 is the handle of something invisible.
         lh = list(ax_12.plot(0, 0)[0:1])
@@ -1316,6 +1350,34 @@ def init(fignum, stats, key0, plot_u, E, P, **kwargs):
         ax_21.set_title("spread. " + sx)
         ax_22.set_title("err. " + sx)
 
+        # Physical grid info (for mapping indices -> physical coords)
+        # arr_mu/arr_xx shape -> (ny, nx)
+        try:
+            ny, nx = arr_mu.shape
+        except Exception:
+            # Fallback: assume square
+            ny = nx = int(np.sqrt(arr_mu.size))
+        Lx, Ly = domainlims
+
+        # Ensure axes limits match physical domain (in units of length)
+        for ax in axs.flatten():
+            ax.set_xlim(0, Lx)
+            ax.set_ylim(0, Ly)
+            # Friendly ticks
+            ax.set_xticks(np.linspace(0, Lx, 5))
+            ax.set_yticks(np.linspace(0, Ly, 5))
+            # relabel in terms of pi if it is periodic
+            if periodic[0]:
+                xticklabels_float = np.linspace(0, domainlims[0] / np.pi, 5)
+                xticklabels_string = xticklabels_float.astype(str)
+                xticklabels_string = np.char.add(xticklabels_string, "π")
+                ax.set_xticklabels(xticklabels_string)
+            if periodic[1]:
+                yticklabels_float = np.linspace(0, domainlims[1] / np.pi, 5)
+                yticklabels_string = yticklabels_float.astype(str)
+                yticklabels_string = np.char.add(yticklabels_string, "π")
+                ax.set_yticklabels(yticklabels_string)
+
         # TODO 7
         # for ax in axs.flatten():
         # Crop boundries (which should be 0, i.e. yield harsh q gradients):
@@ -1370,9 +1432,30 @@ def update(key, E, P):
             #  - ind2sub returns (iy,ix), while plot takes (ix,iy) => reverse.
 
             if ko is not None and not_empty(obs_inds):
-                lh[0] = ax_12.plot(*ind2sub(obs_inds(ko))[::-1], "k.", ms=1, zorder=5)[
-                    0
-                ]
+                inds = obs_inds(ko) if callable(obs_inds) else obs_inds
+                # ind2sub is expected to return (iy, ix) arrays suitable for plotting
+                try:
+                    iy, ix = ind2sub(inds)
+                except Exception:
+                    # If ind2sub returns reversed order, try reversing
+                    tmp = ind2sub(inds)
+                    if len(tmp) >= 2:
+                        iy, ix = tmp[0], tmp[1]
+                    else:
+                        # Give up and plot nothing
+                        iy = ix = np.array([])
+
+                ix = np.asarray(ix)
+                iy = np.asarray(iy)
+
+                # Map integer cell indices to physical cell-centre coordinates
+                if ix.size and iy.size:
+                    x = (ix + 0.5) * (Lx / nx)
+                    y = (iy + 0.5) * (Ly / ny)
+                else:
+                    x = y = []
+
+                lh[0] = ax_12.plot(x, y, "k.", ms=3, zorder=5)[0]
 
             text_t.set_text(format_time(k, ko, t))
 
diff --git a/docs/examples/basic_1.py b/docs/examples/basic_1.py
index c8870794..bd33a713 100644
--- a/docs/examples/basic_1.py
+++ b/docs/examples/basic_1.py
@@ -104,4 +104,4 @@
 # ### Excercise (memory)
 # Why are the replay plots not as smooth as the liveplot?
 #
-# *Hint*: provide the keyword `store_u=True` to `assimilate()` to avoid this.
+# *Hint*: provide the keyword `store_u=True` to `assimilate()` to avoid this.
\ No newline at end of file