2 Data description and image preparation

A thin (ant-farm style) aquarium ( $0.22\, \text{m height}\times 0.445\, \text{m width}\times 0.022\, \text{m depth}$ ) was set up to study the effect of macrofauna on porewater flow. The aquarium was filled with approximately $0.16\, \text{m}$ of near-shore muddy sand type sediment overlain by seawater. After letting the sediment settle for 1 week, a lugworm, Arenicola marina, an abundant macrofauna in the coastal ocean of Europe and North America, was placed in this aquarium. Lugworms are a few centimetres long, and are subsurface deposit feeders who pump water through their burrows into the sediment (e.g. [Reference Riisgård and Banta35]), acting as ecosystem engineers [Reference Volkenborn, Hedtkamp, van Beusekom and Reise47]. Lugworm inject fluid into the sediment by inducing a peristaltic tail-to-head moving pumping wave [Reference Riisgård, Berntsen and Tarp36]. This is done at the feeding pocket, located at the dead end of their J-shaped burrows. When the lugworms are pumping, water is moved through this structure and then injected into the sediment at the dead end. In previous work, A. marina has been modelled with radius of the feeding pocket of 2.5 mm [Reference Meysman, Galaktionov, Gribsholt and Middelburg32]. Given this small size, we describe this as a point source of water originating from outside the modelled domain consistent with observations and brought to the injection location via a burrow that is behind the experimentally observational plane. For the short time interval (on the order of minutes) covered by the images analysed in this study, the burrow structure and injection location are considered constant. Qualitatively, flow velocities decrease with distance form the injection site and depend on the spatial variations in sediment permeability. As lugworms pump intermittently [Reference Woodin, Wethey and Volkenborn54], the resulting flow velocities are temporally variable.

The lugworm was allowed to acclimate for 7 days. During this time, water supply was maintained and air was bubbled in the overlaying water. The aquarium was maintained in the dark.

The flow induced by the lugworm’s activity was visualised by adding fluorescein, a dissolved fluorescent tracer, through a thin pipette to a location in the sediment close to where the burrowing animal resided. Its redistribution in the sediment – caused by pumping of the arenicolid worm – was captured as a series of high-resolution images of the fluorescent flow tracer. The aquaria was illuminated with blue led lights at an approximate wavelength of 450 nm, and a yellow barrier ‘photographic gel’ filter was used to reduce the blue light captured in the images. A 1 g per liter of fluorescein was added to a 32 ppt reagent grade NaCl solution to match the density of the resident pore water (fluid). RGB (red, green, blue) images were captured with a Fujifilm FinePix S5Pro model at resolution $2848\times4256$ pixels every $10\, \text{s}$ . These images are subsequently referred to as real RGB images, and 3 out of the 41 images analysed are displayed in Figure 1.

Figure 1. Real RGB images depicting spatial tracer distribution in a $0.22\, \text{m}\times 0.445\, \text{m}\times 0.022\, \text{m}$ aquarium where bioirrigation is simulated at the relative times $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=190\, \text{s}}$ and $\mathsf{t_2=390\, \text{s}}$ .

The use of a thin aquarium ensures that animal activity leads to a flow visible on the aquarium surface and the assumption of a two-dimensional optical flow. Thus, the incompressible 2D flow is expected to be largely divergence free, except near the location where the fluid is injected, which may not be parallel to the image plane.

The real RGB images are converted to eight-bit grey scale using the linear combination [Reference William51]

(2.1) \begin{equation}\mathsf{F(x,y,t)=0.2989\mathcal{R}(x,y,t)+0.5870\mathcal{G}(x,y,t)+0.1140\mathcal{B}(x,y,t),}\end{equation}

where $\mathsf{\mathcal{R}(x,y,t)}$ , $\mathsf{\mathcal{G}(x,y,t)}$ and $\mathsf{\mathcal{B}(x,y,t)}$ denote red, green and blue colour channel intensities. The coefficients in (2.1) using the Luma rec 601 conversion implemented in the MATLAB rgb2gray command. Arrays arising from (2.1) are subsequently referred to as real grey scale images, and a subset of them are depicted in Figure 2. Importantly, the tracer concentration $\mathsf{C}$ is proportional to $\mathsf{F}$ .

Figure 2. False color gray scale images, whose intensities are obtained using the linear combination in (2.1) with red, green and blue color channels depicted in Figure 1, and cropping boundaries overlaid on top in red.

As we aimed at the quantification of flow in the subsurface, our analyses are carried out on the area delineated in red in Figure 2 only, thus removing the overlying water (top $0.16\, \text{m}$ ) and a small area near the bottom and sides of the aquarium, where the external lighting led to visible variations in brightness that might negatively impact optical flow computations. This image cropping reduces resolution to $1600\times3680$ , resulting in $>11$ million velocity components at the pixel level. To lower computational burden, the real grey scale images were down-sampled by a factor of eight using bilinear interpolation to a resolution of ${200\times 460}$ .

3 Mathematical modeling

Consider the Fokker–Planck equation in (3.1), where $\mathsf{C(x,y,t)}$ denotes tracer concentration, $\mathsf{\textbf{u}(x,y,t)}$ is the velocity vector, and $\mathsf{\nabla=\begin{bmatrix} \mathsf{\frac{\partial}{\partial x}} & \mathsf{\frac{\partial}{\partial y}} \\ \end{bmatrix}^T}$ is the gradient operator [Reference Meysman, Galaktionov, Gribsholt and Middelburg31, Reference Meysman, Galaktionov, Gribsholt and Middelburg32].

(3.1) \begin{equation}\mathsf{(\phi C)_t=\nabla\cdot(\phi D_e \nabla C - \phi{\textbf{u}}C).}\end{equation}

Let $\mathsf{\phi(x,y,t)}$ denote sediment porosity, the ratio of void to total volume, on the order of $0.7$ in sandy sediments inhabited by Arenicola marina. Here, it is assumed to be a spatio-temporally constant (for simplicity), so that

(3.2) \begin{equation}\mathsf{C_t=\nabla\cdot\left(D_e \nabla C - {\textbf{u}}C\right).}\end{equation}

The scalar $\mathsf{D_e\in\mathbb{R}}$ in (3.2) denotes the effective diffusion of the sediment, and is given by

(3.3) \begin{equation} \mathsf{D_e=\frac{D}{1-2\ln\!{\phi}}},\end{equation}

where $\mathsf{D\in\mathbb{R}}$ is the molecular diffusion coefficient in aqueous solution. For fluorescein, $\mathsf{D=4.25\times 10^{-10}\, \text{m}^2\, \text{s}^{-1}}$ so that $\mathsf{D_e=2.48\times 10^{-10}\, \text{m}^2\, \text{s}^{-1}}$ [Reference Cullbertson, Jacobson and Michael Ramsey13]. Consider the Péclet number given by

(3.4) \begin{equation}\mathsf{\text{Pe}=\frac{\ell\|{\textbf{u}}\|_2}{D_e},}\end{equation}

where $\mathsf{\ell>0}$ is a length scale and $\mathsf{\|\cdot\|_2}$ denotes a Euclidean norm. Advection is the more prominent solute transport mechanism when $\mathsf{\text{Pe}\gg 1}$ . Now $\ell$ is $\mathsf{\mathcal{O}\!\left(10^{-2}\right)\, \text{m}}$ , in line with the dimensions of the aquarium from which the real RGB images arise. Also, previous research suggests that the flow speed $\mathsf{\|{\textbf{u}}\|_2}$ of water induced by bioirrigating Arenicola marina is approximately $\mathsf{\mathcal{O}\big(10^{-6}\big)\, \text{m}\, \text{s}^{-1}-}$ $\mathsf{\mathcal{O}\big(10^{-5}\big)}\, \text{m}\, \text{s}^{-1}$ [Reference Kaza22, Reference Volkenborn, Polerecky, Wethey and Woodin48]. This, alongside the earlier definitions of porosity and diffusion coefficient, indicates that $\mathsf{\text{Pe}\gg 1}$ . This motivates the assumption that diffusion is negligible, particularly in the region of interest in the vicinity of the burrow, so that (3.2) can be rewritten as

(3.5) \begin{equation} \mathsf{-C_t=\nabla\cdot({\textbf{u}}C).}\end{equation}

We will assume incompressibility and the condition $\mathsf{\nabla\cdot{\textbf{u}}}=0$ , whose validity for the current context was discussed in previous research [Reference Kaza22, Reference Volkenborn, Polerecky, Wethey and Woodin48]. It then follows that

(3.6) \begin{equation}\begin{aligned} \mathsf{-C_t}&=\mathsf{C \nabla\cdot{\textbf{u}}+{\textbf{u}}\cdot\nabla C.}\\ &=\mathsf{{\textbf{u}}\cdot\nabla C.}\end{aligned}\end{equation}

The connection between (4.1) and (3.6) are made by assuming that tracer concentration is proportional to grey scale intensity (i.e. $\mathsf{C}\propto\mathsf{F}$ ) [Reference Kaza22].

4 Horn–Schunck method

For a velocity field to be estimated, a mathematical expression relating it to grey scale intensity is required. Consider a pair of images, whose grey scale intensities are defined by $\mathsf{F_0(x,y,t_0)}$ and $\mathsf{F_1(x,y,t_1)}$ for $\mathsf{t_0<t_1}$ . An optical flow $\mathsf{\textbf u}(x,y,t)=[\mathsf{u(x,y,t)}\enspace\mathsf{v(x,y,t)}]^T$ is an apparent motion field mapping one image to the next. Conventional optical flow problem relies on the (pixel) brightness constancy approximation. This is expressed mathematically in (4.1), where the subscripts $\mathsf{x}$ , $\mathsf{y}$ and $\mathsf{t}$ denote derivatives with respect to the $\mathsf{x}$ -axis, $\mathsf{y}$ -axis and time [Reference Horn and Schunck21].

(4.1) \begin{equation}\mathsf{F_t+uF_x+vF_y=0}.\end{equation}

In the current context, this approximation can be justified by requiring the null divergence of the velocity in (3.5) at every image pixel. This is, however, ill-posed since two velocity components have to be determined from one equation. To address this issue, a classical Horn–Schunk method was developed based on a variational principle [Reference Horn and Schunck21]. In particular, the Horn–Schunck method further constrains the optical flow to be spatially smooth, which can be expressed mathematically as

(4.2) \begin{equation} \mathsf{\|\nabla u\|_2^2+\|\nabla v\|_2^2<\infty.}\end{equation}

Locating the optical flow that minimises (4.1) subject to (4.2), with respect to the spatial region $\mathsf{\Omega}$ that the images depict, can be expressed mathematically as

(4.3) \begin{equation} \mathsf{\textbf{u}}=\underset{{\textbf{u}}}{\text{argmin}}\int_{\mathsf{\Omega}} \underbrace{\left(\mathsf{\textbf{u}}\cdot\nabla \mathsf{F}+\mathsf{F}_{\mathsf{t}}\right)^{\mathsf{2}}+\mathsf{\lambda}^{\mathsf{2}}\left(\|\nabla \mathsf{u}\|_{\mathsf{2}}^{\mathsf{2}}+\|\nabla \mathsf{v}\|_{\mathsf{2}}^{\mathsf{2}}\right)}_{\mathsf{\mathcal{I}}} \mathsf{d\Omega}\end{equation}

where $\lambda\in\mathbb{R}$ is a smoothing parameter. Following [Reference Horn and Schunck21], the Euler–Lagrange equations for (4.3) are given by the following coupled system of partial differential equations (PDEs), where $\mathsf{\Delta}$ is the Laplacian operator,

(4.4) \begin{equation}\begin{aligned} \mathsf{F_x\!\left(F_xu+F_yv+F_t\right)} & \,\mathsf{=}\, \mathsf{\lambda^2\Delta u,}\\ \mathsf{F_y\!\left(F_xu+F_yv+F_t\right)} & \,\mathsf{=}\, \mathsf{\lambda^2\Delta v.}\end{aligned}\end{equation}

At the image pixel level, the system in (4.4) can be numerically solved in which the Laplacian is approximated by

(4.5) \begin{equation} \mathsf{\Delta u\approx \overline{\overline {u}}-u},\qquad\overline{\overline {u}}=\frac{1}{12}\begin{bmatrix}1 & \quad 2 & \quad 1\\[2pt]2 &\quad 0 & \quad 2\\1 & \quad 2 & \quad 1\\\end{bmatrix}\ast u.\end{equation}

The $\ast$ in (4.5) means a 2D convolution on the image pixel. In other words, $\overline{\overline {u}}$ is a local-weighted average. After substituting the Laplacian approximation in (4.5) into (4.4), we can apply the Jacobi iteration to obtain the following iterative scheme, where the subscript $\mathsf{i}$ denotes iteration indices,

(4.6) \begin{equation}\begin{aligned} &\mathsf{u_{i+1}=\overline{\overline{u}}_i-\frac{F_x\!\left(F_x\overline{\overline{u}}_i+F_y\overline{\overline{v}}_i+F_t\right)}{\lambda^2+F_x^2+F_y^2}.}\\[3pt] &\mathsf{v_{i+1}=\overline{\overline{v}}_i-\frac{F_y\!\left(F_x\overline{\overline{u}}_i+F_y\overline{\overline{v}}_i+F_t\right)}{\lambda^2+F_x^2+F_y^2}.}\end{aligned}\end{equation}

The differential nature of the Horn–Schunk method can make it difficult to estimate a large displacement and not robust to noise. In practice, a multi-resolution frame work (e.g. see [Reference Black and Anandan5, Reference Ruhnau, Kohlberger, Schnörr and Nobach37]) is commonly used to mitigate these issues. This approach is based on an incremental estimate that refines a previous estimate using finer resolution image at the next higher resolution level, which is commonly referred to as ‘pyramid’ level. At the top of the pyramid level, the estimation is carried out at the original resolution. At the lower levels, the image is repeatedly down-sampled for a pre-specified factor. The details of the multi-resolution Horn–Schunk method is discussed in Appendix A. The number of pyramid levels is chosen as $\mathsf{\mathcal{L}=3}$ here to obtain the results below. Ideally, magnitudes of additional terms added to a current correction increment in (A3) should depend on how faithfully the associated optical flow mimics the unknown true solution. So if a current correction increment produces an optical flow that almost exactly mimics the true solution, then magnitudes of additional terms in (A3) at subsequent iterations should be approximately zero. The magnitudes of the quotients in (A4) can be controlled by the value of $\mathsf{\lambda}$ .

The smoothing parameter should not be so ‘small’ that these quotients are considered ‘big’ in regions where contrast is low. Otherwise, optical flow speeds could then temporally increase to a point where they are viewed as being unrealistically ‘large’. However, the smoothing parameter should not be too ‘large’. This could induce an unrealistic high level of spatial smoothing. By trial and error, the value of $\mathsf{\lambda=50}$ is chosen because it appears to be one of the ‘smaller’ smoothing parameter values where flow speeds do not become unrealistically ‘large’ over time. This is examined by plotting a frequency distribution of flow speeds at each relative time point in Figure 3.

Figure 3. Each column shows the frequency distribution of optical flow speeds obtained by the multi-resolution Horn–Schunck method for $\mathsf{\lambda=50}$ and $\mathsf{M=1}$ .

The number of iterations of (A3) at each pyramid level, denoted by $\mathsf{M}$ , also needs to be defined. The results here are shown for $\mathsf{M=1}$ as per the following justification. It appears that when using $\mathsf{M>1}$ , at later relative time points, the only region where flow points downwards is between the likely location of the organism and the bottom edge of the down-sampled and cropped real grey scale images. This suggests that there is outflow and no inflow at the top boundary of these images, which roughly coincides with the sediment–water interface according to Figure 2. However, given that the sides and bottom of the aquarium are impermeable, net flow across the sediment-water interface should be zero. We expect some outflow in the regions near the aquarium walls were cropped due to uneven lighting (Section 2). However, the qualitative assessment of tracer distribution in the overlying water (Figure 2), indicates that this flow is limited. Thus, $\mathsf{M}$ is chosen to enforce a little inflow at later time points.

Figure 4 reveals undesirably ‘large’ divergence magnitudes where they ought to be comparatively ‘small’, i.e., away from the location where the worm activity induces flow (see spreading arrows near the center of the top left panel). Furthermore, largest velocities are computed away from this point, with no clear source of fluid flow being discernible. Thus, additional ways of imposing that divergence magnitudes immediately tend to zero upon departure from the point of injection are investigated. This is developed in the next section.

Figure 4. The first column depicts optical flows, denoted by white arrows, arising from the multi-resolution Horn–Schunck method when $\mathsf{\lambda=50}$ and $\mathsf{M=1}$ at the relative times $\mathsf{t_0=0\, \text{s}}$ (first row), $\mathsf{t_1=190\, \text{s}}$ (second row) and $\mathsf{t_2=390\, \text{s}}$ (third row). They are overlaid on top of the down-sampled and cropped real grey scale images associated with those relative times. The second column depicts the associated divergences. Note that all parts of this figure have been further down-sampled by a factor of four, for ease of presentation.

5 Derivation of parametric framework

Arenicola marina inject water into the sediment at the dead end of often J-shaped burrows. Simulations of bioirrigation by the lugworm have shown that flow fields are largely unaffected by the burrow structure itself [Reference Meysman, Galaktionov, Gribsholt and Middelburg31]. Hence, we consider flow in the sediment being driven by a fluid source located at a small injection pocket. Away from injection pocket, flow is considered two-dimensional, parallel to the imaging surface, due to closely spaced front and back wall of the thin aquarium. Therefore, the condition $\mathsf{\nabla\cdot{\textbf{u}}=0}$ in incorporated within the Bayesian framework through the prior distribution with the exception to this at the point $\mathsf{(a,b)}$ , where the organism injects fluid. In particular, we will write an ansatz for the velocity potential $\mathsf{\Phi(x,y,t)}$ of a two-dimensional incompressible point source with strength $\mathsf{c}>0$ centred on the point $\mathsf{(a,b)\in\mathbb{R}^2}$ , subject to boundary conditions as discussed in Appendix A. The injection point induces fluid expansion with strength $\mathsf{c}$ , suggesting that $\mathsf{\nabla\cdot{\textbf{u}}=c}$ . All above requirements suggest the condition

(5.1) \begin{equation} \mathsf{\nabla\cdot{\textbf{u}}=c\delta(x-a,y-b).}\end{equation}

We can write (5.1) in terms of a velocity potential by using $\mathsf{{\textbf{u}}=\nabla \Phi}$ .

(5.2) \begin{equation} \mathsf{\Delta \Phi=c\delta(x-a,y-b).}\end{equation}

To deal with the delta function on the right-hand side of (5.2), we follow the method demonstrated in [Reference Batchelor2] to approximate the solution of (5.2) with that of

(5.3) \begin{equation}\mathsf{\Delta \Phi=\nu\Phi},\end{equation}

where $\mathsf{\nu\in\mathbb{R}}$ is an eigenvalue. The ansatz of velocity potential of a two-dimensional incompressible point source with strength $\mathsf{c}$ centred on the location $\mathsf{(a,b)}$ is given by

(5.4) \begin{equation}\mathsf{\Phi(x,y)=\frac{c}{2\pi}\ln{\sqrt{\left((x-a)^2+(y-b)^2 \right)}}.}\end{equation}

We use the separation of variables to solve (5.4), see Appendix B for the outline of the derivation,

(5.5) \begin{equation}\begin{aligned}\mathsf{\Phi(x,y)}&\,\mathsf{=-\sum_{j=0}^\infty\sum_{k=0}^\infty\frac{4c}{L_xL_y\!\left[ \left(\frac{k\pi}{L_x}\right)^2+\left(\frac{(1+2j)\pi}{2L_y}\right)^2 \right]}\cos\!{\left( \frac{k\pi}{L_x}\, a \right)}\cos\!{\left(\frac{(1+2j)\pi}{2L_y}\, b \right)}}\\&\mathsf{\, \, \, \times\cos\!{\left( \frac{k\pi}{L_x}\, x \right)}\cos\!{\left(\frac{(1+2j)\pi}{2L_y}\, y \right)},}\end{aligned}\end{equation}

where $\mathsf{[0,L_x]\times[0,L_y]}$ is the area under computation, see Figure B.1. The infinite sums in (5.5) is approximated with the finite sums with the first $\mathsf{M_x}$ terms for the inner sum and $\mathsf{M_y}$ for the outer sum. In what follows, we choose $\mathsf{M_x=M_y=200}$ . Utilising the ansatz in (5.5) means that once the parameters $\mathsf{\{a,b,c\}}$ are estimated, flow at each grid cell can be obtained using the relation $\mathsf{{\textbf{u}}=\nabla\Phi}$ .

6 MCMC for flow field estimation

The relation $\mathsf{{\textbf{u}}=\nabla\Phi}$ can be used to rewrite (4.1) in terms of velocity potentials [Reference Basnayake1, Reference Horn and Schunck21, Reference Luttman, Bollt, Basnayake, Kramer and Tufillaro25].

(6.1) \begin{equation} \mathsf{-F_t=F_x\Phi_x+F_y\Phi_y.}\end{equation}

Substituting the finite approximation of (5.5) into (6.1) yields

(6.2) \begin{align}\mathsf{-F_t}&\mathsf{\approx\sum_{j=0}^{M_y}\sum_{k=0}^{M_x} \frac{4c}{L_xL_y\nu^{jk}}\cos\!{\left( \frac{k\pi}{L_x}\, a \right)}\cos\!{\left(\frac{(1+2j)\pi}{2L_y}\, b \right)}\bigg[{-}\frac{F_xk\pi}{L_x}\sin\!{\left( \frac{k\pi}{L_x}\, x \right)}\cos\!{\left(\frac{(1+2j)\pi}{2L_y}\, y \right)}}\nonumber\\&\mathsf{\, \, \, -\frac{F_y(1+2j)\pi }{2L_y}\cos\!{\left( \frac{k\pi}{L_x}\, x \right)}\sin\!{\left(\frac{(1+2j)\pi}{2L_y}\, y \right)}\bigg].}\end{align}

This acts as an observation system relating the parameters $\mathsf{\{a,b,c\}}$ contained within the state vector to observed negated temporal grey scale intensity derivatives. Let $\mathsf{\textbf{z}:=-F_t}$ and $\mathsf{\textbf{x}:=[a,b,c]^T}$ and abbreviate (6.2) by

(6.3) \begin{equation}\mathsf{\textbf{z}=\mathcal{H}(\textbf{x})}.\end{equation}

We will investigate the uncertainty of the parameters for the three observations of $\mathsf{\textbf{z}:=-F_t}$ approximated from the (grey-scale) image pairs at the relative time point $(t=0,t=10)$ , $(t=10,t=20),\ldots,(t=40,t=50)$ s. These will be denoted by $\mathsf{\textbf{z}(0),\textbf{z}(10),\ldots,\textbf{z}(40)}$ , respectively. We will assimilate the observation one at a time to sequentially sampling $\mathsf{p(\textbf{x}|\textbf{z(0)})}$ , $\mathsf{p(\textbf{x}|\textbf{z(0)},\textbf{z(10)})}$ ,…, $\mathsf{p(\textbf{x}|\textbf{z(0)},\ldots,\textbf{z(40)})}$ .

6.1 Setting likelihood function and prior distribution

To adopt a Bayesian framework, we need to assume a conditional distribution of $\mathsf{\textbf{z}}$ given $\mathsf{\textbf{x}}$ and the prior distribution of $\mathsf{\textbf{x}}$ . For simplicity, we assume the likelihood function of the observation $\textbf{z}$ at each time point to have the kernel (i.e. the non-normalised part of the distribution)

(6.4) \begin{equation}\mathsf{p({\textbf{z}}|{\textbf{x}})\propto\exp\!{\left(-\frac{1}{2}\!\left({\textbf{z}}-\mathcal{H}({\textbf{x}})\right)^TR^{-1}\!\left({\textbf{z}}-\mathcal{H}({\textbf{x}})\right)\right)},}\end{equation}

where $\mathsf{R=\sigma_2^2 I}$ and $\mathsf{I}$ is the identity matrix of an appropriate size. An empirical value of $\mathsf{\sigma_2^2=1\times10^{-3}\, \text{s}^{-2}}$ is adopted here following a numerical experiment on synthetic data with similar features to our current real-world data [Reference McCormack28]. We use the same likelihood function for $\mathsf{\textbf{z}}$ at all time points. The state vector now contains source strengths, and $\mathsf{x}$ - and $\mathsf{y}$ -coordinates, which have different physical properties. For the source $\mathsf{x}$ -coordinates, the horizontal distance from the source to the right-hand boundary $\mathsf{\widetilde a\in\mathbb{R}}$ is to be estimated instead of $\mathsf{a}$ , where $\mathsf{\widetilde a=L_x-a}$ . We treat these parameters as being probabilistically independent, i.e.,

(6.5) \begin{equation} \mathsf{p({\textbf{x}})=p(\widetilde a)p(b)p(c).}\end{equation}

The assumption of independence is justified because the organism’s location in the sediment should not affect the flow it induces, nor does the burrowing depth impact the lateral positioning.

We now describe the choice of the prior distribution at the first time point and will subsequently discuss how the prior for later time points are constructed. The marginal prior distribution for the source strength is assumed to be a Gamma distribution to ensure the positivity. The prior assumption $\mathsf{E[c]=2\times 10^{-5}\, \text{m}^2\, \text{s}^{-1}}$ is made because flow speeds residing approximately within the interval $\mathsf{\mathcal{O}\big(10^{-6}\big)\, \text{m}\, \text{s}^{-1}-}$ $\mathsf{\mathcal{O}\big(10^{-5}\big)\, \text{m}\, \text{s}^{-1}}$ , which appear to be induced by bioirrigating Arenicola marina in previous research [Reference Kaza22, Reference Volkenborn, Polerecky, Wethey and Woodin48]. The marginal prior variance is taken to be $\mathsf{\text{Var}[c]=2.5\times 10^{-11}\, \text{m}^4\, \text{s}^{-2}}$ . It should be ‘large’ enough that it can be reduced when observations are assimilated, whilst not being so ‘large’ that the information quantified by the marginal prior mean is overpowered by observations. The marginal prior means are derived from the $\mathsf{x}$ -and $\mathsf{y}$ -coordinates of the pixel in the first down-sampled and cropped real grey scale image with the greatest intensity, because intuitively this seems to be a realistic location of fluid injection. Let $\mathsf{\delta_x=L_x^{-1}}$ and $\mathsf{\delta_y=L_y^{-1}}$ . The marginal prior distributions at the first relative time point for the parameters $\mathsf{b}$ is assumed to be a Gaussian distribution with mean $\mathsf{E[b]=0.07\, \text{m}}$ and variance $\mathsf{\text{Var}[b]=\frac{5\delta y^2}{2}\, \text{m}^2}$ . Finally, consider the parameter $\mathsf{\widetilde a}$ . A negative value of $\mathsf{\widetilde a}$ is unrealistic. As such, a Gamma marginal prior distribution is employed for $\mathsf{\widetilde a}$ . So we assume a Gamma distribution with mean $\mathsf{E[\widetilde{a}]=0.2058\, \text{m}}$ and variance $\mathsf{9\delta x^2\, \text{m}^2}$ .

Note that the prior mean is set to be ‘small’, with a view to quantifying the fact that the source lies very close to the right-hand boundary. Their related standard deviations of $\mathsf{p(\widetilde a)}$ and $\mathsf{p(b)}$ are ‘small’ compared with the width and height of the aquarium to ensure that the source $\mathsf{x}$ -and $\mathsf{y}$ -coordinates recovered from the Metropolis–Hastings (MH) sampler remain ‘close’ to the regions with the greatest intensity. With above prior assumption, the MH sampler can be used to sample the posterior distribution

(6.6) \begin{equation} \mathsf{p(\textbf{x}|\textbf{z(0)})\propto p(\textbf{z(0)}|\textbf{x})p(\widetilde a)p(b)p(c)}.\end{equation}

The posterior distribution of the parameter vector $\mathsf{\textbf{x}}$ given $\mathsf{\textbf{z(0)}}$ could then be used as the prior density for the next observation $\mathsf{\textbf{z(10)}}$ . However, we only have a sample of the posterior distribution, not the analytical form. To deal with this issue, we assume the independent prior for $\mathsf{\textbf{x}}$ when assimilating $\mathsf{\textbf{z(0)}}$ and then construct the prior distribution as above by using the posterior (marginal) mean and variance for each parameter $\mathsf{a,b,c}$ given $\mathsf{\textbf{z(0)}}$ instead. The process is then repeated for the rest of the observations.

6.2 Metropolis–Hasting

We use the MH sampler [Reference Hastings19] to sample (6.6). The MH algorithm involves accepting or rejecting a proposed sample with the acceptance probability

(6.7) \begin{equation} \mathsf{\gamma=\min\!\left\{1,\frac{p\!\left({\textbf{x}}^\ast|{\textbf{z}}\right)}{p\!\left({\textbf{x}}_i|{\textbf{z}}\right)}\frac{p\!\left({\textbf{x}}_i|{\textbf{x}}^\ast\right)}{p\!\left({\textbf{x}}^\ast|{\textbf{x}}_i\right)}\right\},}\end{equation}

where $\textbf{x}^\ast$ is a (new) proposed sample and $p(\textbf{x}^\ast|\textbf{x}_i)$ is a proposal distribution of $\textbf{x}^\ast$ given the current sample $\textbf{x}_i$ . We assume the proposal distribution of the following form

(6.8) \begin{equation}\begin{aligned}\mathsf{p\!\left({\textbf{x}}^{\ast}|{\textbf{x}}_{i}\right)}&\, \mathsf{=p\!\left(\widetilde a^{\ast},b^{\ast},c^{\ast}|\widetilde a_{i},b_{i},c_{i}\right)}\\&\mathsf{=p\!\left(\widetilde a^{\ast}|\widetilde a_{i}\right)\!p\!\left(b^{\ast}|b_{i}\right)\!p\!\left(c^{\ast}|c_{i}\right).}\end{aligned}\end{equation}

The marginal proposal distribution $\mathsf{p\!\left(\widetilde a^{\ast}|\widetilde a_{i}\right)}$ is chosen to be the Gamma distribution with mean $\widetilde a_{i}$ and variance is obtained by pre-multiplying marginal prior variance of $\widetilde a$ by a factor of 0.36 here. The marginal proposal distribution $\mathsf{p\!\left(\widetilde c^{\ast}|\widetilde c_{i}\right)}$ is also defined by the Gamma distribution with mean $\widetilde c_{i}$ and variance is obtained as above. A Gaussian marginal proposal distribution is defined for source $\mathsf{y}$ -coordinates to enforce that injection can realistically take place either above or below the current chain location. In particular, $\mathsf{p(b^{\ast}|b_{i})}$ is a Gaussian distribution with mean $b_i$ and variance is the marginal prior variance of b multiplied by 0.36. The ad hoc choice of the above pre-multiplying factor is chosen from a trial and error to strike a balance between the acceptance and convergence rates; when the proposal variance is too ‘small’, only ‘small’ steps can be taken across the state space and convergence is slow. However, we make no attempt to optimally tune this factor here.

6.3 Convergence diagnosis

In what follows we will run multiple Markov chains, starting in different regions of the state space, and consider both inter-chain mixing and intra-chain mixing for convergence diagnosis. Suppose that there are $\mathsf{\mathcal{C}\in\mathbb{N}}$ Markov chains. Each chain ends up of length $\mathsf{N_{spin}+N}$ , where the first $\mathsf{N_{spin}}$ particles constitute an initial transient called a spin-up period. The latter $\mathsf{N}$ particles constitute the ensemble that the Metropolis–Hastings sampler generates. The variance of the $\mathsf{j}$ -th Markov chain (for $\mathsf{j\in\{1,\dots,\mathcal{C}\}}$ ) is given by

(6.9) \begin{equation}\mathsf{s_j^2=\frac{1}{\widetilde{N}-1}\sum_{i=1}^{\widetilde{N}}\!\left(c_{ij}-\overline c_j\right)^2}\end{equation}

after the $\mathsf{\widetilde{N}}$ -th particle is obtained. The terms $\mathsf{c_{ij}}$ and $\mathsf{\overline c_j}$ denote the $\mathsf{i}$ -th member of the $\mathsf{j}$ -th Markov chain $\big($ for $\mathsf{i\in\big\{1,\dots,\widetilde{N}\big\}}\big)$ , and the mean of the $\mathsf{j}$ -th Markov chain (respectively). Note that $\mathsf{s_j^2}$ is written in terms of source strengths but it can be computed for source $\mathsf{x}$ - and $\mathsf{y}$ -coordinates as well. The expression in (6.9) is used to define the within-chains variance, which is given by [Reference Gelman and Rubin18]

(6.10) \begin{equation}\mathsf{\mathcal{W}=\frac{1}{\mathcal{C}}\sum_{j=1}^{\mathcal{C}}s_j^2.}\end{equation}

The between-chains variance, denoted by $\mathcal{B}$ , is defined by [Reference Gelman and Rubin18, Reference Lambert23]

(6.11) \begin{equation} \mathsf{\frac{\mathcal{B}}{{\widetilde{N}}}=\frac{1}{\mathcal{C}-1}\sum_{j=1}^\mathcal{C}\!\left( \overline{c}_j-\overline{c}\right)^2.}\end{equation}

The term $\mathsf{\overline{c}}$ in (6.11) denotes the overall mean across all $\mathsf{\mathcal{C}}$ Markov chains of source strengths. As suggested in [Reference Gelman and Rubin18, Reference Lambert23], the following ratio

(6.12) \begin{equation} \mathsf{\widehat{R}=\sqrt{\frac{\mbox{Var}[\widehat c|{\textbf{z}}]}{\mathcal{W}}}}\end{equation}

could be used to judge when a set of Markov chains converges, where

(6.13) \begin{equation}\mathsf{\mbox{Var}[\widehat c|{\textbf{z}}]=\left(1-\frac{1}{\widetilde{N}}\right)\!\mathcal{W}+\frac{\mathcal{B}}{\widetilde{N}}.}\end{equation}

Typically, it is suggested that the MCMC should be run until $\mathsf{\widehat{R}<1.1}$ [Reference Lambert23].

6.4 Results: MCMC

To avoid cluttering the presentation of the result, results are shown only at the relative times $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=20\, \text{s}}$ and $\mathsf{t_2=40\, \text{s}}$ . This experiment runs 6 Markov chains and the initial positions of these chains are independently drawn from the prior distribution explained in Section 6.1. We decide on using six chains after the fact that we obtained an acceptable value of $\mathsf{\widehat{R}}$ in this experiment. The spin-up period consists of the first 500 particles and the overall length of each chain, including the spin-up period, is 2000.

Figure 5 depicts evolving Markov chain means for source strengths, $\mathsf{x}$ -coordinates and $\mathsf{y}$ -coordinates. For each panel, chain means appear to converge to similar regions of the state space. Figure 6 depicts Markov chain distributions of source strengths, as well as its $\mathsf{x}$ - and $\mathsf{y}$ -coordinates, at the final iteration (excluding the spin-up period) for the relative times points $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=20\, \text{s}}$ and $\mathsf{t_2=40\, \text{s}}$ . Their shapes, for each parameter–time combination, appear to be roughly the same despite some noise. However, this is not always an evidence of ‘good’ intra-chain mixing because the marginal prior standard deviations defined earlier for $\mathsf{b}$ and $\mathsf{\widetilde{a}}$ appear to be ‘small’ compared with the width and height of the aquarium under consideration. For a particular parameter–time combination, one could argue that Markov chain distributions mimic each other due to their starting positions being too close together instead of because they all remain invariant when more particles are obtained.

Figure 5. Evolving Markov chain means for source strengths (first column), $\mathsf{x}$ -coordinates (second column) and $\mathsf{y}$ -coordinates (third column). The relative times $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=20\, \text{s}}$ and $\mathsf{t_2=40\, \text{s}}$ constitute the first, second and third rows (respectively). The dashed lines denote the end of the spin-up period, which constitutes the first 500 particles.

Figure 6. Markov chain distributions at the final iteration for source strengths (first column), $\mathsf{x}$ -coordinates (second column) and $\mathsf{y}$ -coordinates (third column). The relative times $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=20\, \text{s}}$ and $\mathsf{t_2=40\, \text{s}}$ constitute the first, second and third rows (respectively). Note that the particles constituting the spin-up period are ignored.

The within- and between-chains variances are now examined. The evolving $\mathsf{\widehat{R}}$ ratios are presented in Figure 7 at the relative times $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=20\, \text{s}}$ and $\mathsf{t_2=40\, \text{s}}$ . These appear to satisfy the condition $\mathsf{\widehat{R}<1.1}$ as suggested in [Reference Lambert23] upon reaching the final iteration. This, alongside Figures 5 and 6, reasonably suggest that Markov chains have mixed and converged to their relevant marginal posterior distributions.

Figure 7. Evolving $\mathsf{\widehat{R}}$ ratios for source strengths (first column), $\mathsf{x}$ -coordinates (second column), and $\mathsf{y}$ -coordinates (third column) at the relative times $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=20\, \text{s}}$ and $\mathsf{t_2=40\, \text{s}}$ . The vertical dashed line denotes the end of the spin-up period. The horizontal dashed lines denote the threshold $\mathsf{\widehat{R}<1.1}$ suggested in [Reference Lambert23].

There appear to be a number of situations in Figures 8 and 9 where within- and between-chains variances stabilise after the end of the spin-up period. There seem to be some exceptions to this as well. For example, the between-chains variances of source strengths at the relative time $\mathsf{t_2=40\, \text{s}}$ rapidly increase at later iterations. One possible explanation for this behaviour in the bottom left-hand panel of Figure 9 is that the posterior variance reduction at earlier times could potentially produce ‘small’ marginal prior variances. Starting positions of Markov chains could end up being so ‘close’ together that their evolving means begin to diverge slightly after the spin-up period has elapsed.

Figure 8. Evolving within-chains variances in the log–log scale for source strengths (first column), $\mathsf{x}$ -coordinates (second column), and $\mathsf{y}$ -coordinates (third column) at the relative times $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=20\, \text{s}}$ , and $\mathsf{t_2=40\, \text{s}}$ . The dashed line denotes the end of the spin-up period.

Figure 9. Evolving between-chains variances in the log–log scale for source strengths (first column), $\mathsf{x}$ -coordinates (second column) and $\mathsf{y}$ -coordinates (third column) at the relative times $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=20\, \text{s}}$ and $\mathsf{t_2=40\, \text{s}}$ . The dashed line denotes the end of the spin-up period.

To examine temporal dynamics of optical flows, the flow speed at every grid cell is computed. The frequency distributions of flow speeds at each relative time point in Figure 10 are then constructed and compared with those computed by the multi-resolution Horn–Schunk. These flow speeds are consistent with previous research [Reference Kaza22, Reference Volkenborn, Polerecky, Wethey and Woodin48].

Figure 10. (Left) The frequency distributions of optical flow speeds at the first 5 image pairs computed using the multi-resolution Horn–Schunk, see again Figure 3. (Right) The frequency distributions of optical flow speeds computed using the finite approximation of the ansatz in (5.5), the relation $\mathsf{{\textbf{u}}=\nabla\Phi}$ and the parameter estimates in Table 1.

The results from the HS and our MCMC approach differ substantially, with a temporally more variable and slower flow being computed using HS (Figure 10). Active bioirrigation by Arenicola marina typically occurs over periods of 15 min or longer (e.g. [Reference Timmerman, Banta and Glud46, Reference Wohlgemuth, Taylor and Grieshaber52]). This suggests that the estimates of our MCMC approach are superior to those obtained using HS, where pumping is estimated to vary substantially over a period of less than a minute. Yet porewater pressure time series measurements reveal variations at high frequencies during periods of pumping, and have revealed strong variations associated with burrowing and backward pumping [Reference Volkenborn, Polerecky, Wethey and Woodin48]. Thus, while less common and hence less likely to be captured in our image data, variations in flow magnitudes on timescales shorter than minutes as indicated by the HS approach cannot be ruled out.

Table 1. Markov chain means and standard deviations with respect to both particles and chains (excluding the spin-up period) for source strengths, and its $\mathsf{x}$ -and $\mathsf{y}$ -coordinates. These are obtained using the Metropolis–Hastings sampler and the first five pairs of consecutive down-sampled and cropped real grey scale images

Contrary to what is expected based on our understanding of the biology of Arenicola marina which injects fluid at its dead-end burrow, the HS approach shows fastest flow located further and further away from the location of the burrowing organism over time (Figure 11), largely tracking the gradients in the tracer concentration field. In these distal gradient regions, flow is well constrained by the observational data. The source of flow is the pumping by Arenicola marina, situated near the injection site. Such an injection site is reflected in the location of the divergence in the MCMC approach, but not in the widely distributed divergence in the HS approach. Finally, with a localised injection site, flow velocities should decrease away from the source as the fluid is spreading out. The MCMC approach correctly simulates these properties, whereas the HS method does not. Thus, our work highlights the importance of the data analysis method when inferring the flow field. In this example, the comparison with biological knowledge clearly indicates the advantage of our MCMC approach over traditional methods.

Figure 11. Comparison of the optical flow and divergence between using the finite approximation of the ansatz in (5.5) and Horn–Schunk method in Section 4. The plot shows only the cropped area where Arenicola marina is active. The white arrows in the first and second columns denote optical flows. From top to bottom row, the results are shown at the relative times $\mathsf{t_0=0\, \text{s}}$ , $\mathsf{t_1=20\, \text{s}}$ and $\mathsf{t_2=40\, \text{s}}$ , respectively. They are overlaid on top of the down-sampled and cropped real grey scale images associated with those relative times. The third and fourth columns depict the associated divergences. Note that the first and second columns share the same colour scale and likewise for the third and fourth columns.

7 Summary and outlook

This work develops a parametric model to estimate the flow velocity in the context of bioirrigation. A localised divergence-free flow velocity is desired in the vicinity of the injection point, which is difficult to achieve with the conventional flow velocity framework method such as the Horn–Schunk method. We demonstrated that the issue of undesired divergence of the flow velocity can be significantly mitigated by using the proposed model in combination with MCMC. The efficacy of the MCMC method in this work indicates the superiority compared to other methods when applied to general problems associated with localised sources of flow/pressure field generation.

Despite the improved result in this initial work, several challenges remain open to overcome. In particular, the model (3.5) can be inaccurate further away from the injection site where the effect of diffusion is considerable and has to be taken into account in the model. This must be addressed if flow velocity estimation has to be extended to the entire area of the image presented in this work. Taking into account diffusion, or considering reactive tracers that inform our understanding of the biogeochemical turnover and efflux of, for example, greenhouse gases such as $CO_2$ , $N_2O$ or $CH_4$ from the subsurface would add more parameters, rendering MCMC infeasible in practice. This should motivate a more advanced (approximate) inference to deal with this inverse problem, e.g., iterative Ensemble Kalman filter (IEnKF) for solving inverse problems [Reference Chada, Iglesias, Roininen and Stuart8, Reference Emerick and Reynolds15, Reference Sakov, Haussaire and Bocquet38]. Another direction to further improve the estimate by setting up a model to allow a temporal-smoothing parameter into the framework instead of using the posterior distribution as a prior distribution in a subsequent image pair. This will in effect avoid ‘cutting the feedback’ issue discussed in [Reference Plummer33].