Global existence, blow-up behavior and numerical simulations for a class of chemotaxis-driven fish-mussel systems

Technical Summary: Global Existence, Blow-Up Behavior, and Numerical Simulations for a Class of Chemotaxis-Driven Fish-Mussel Systems

Problem Statement
This work investigates a chemotaxis-driven ecosystem model describing the spatiotemporal dynamics of nutrients ($u$), fish populations ($v$), and mussel populations ($w$) within a bounded domain $\Omega \subset \mathbb{R}^n$. The model is formulated as a coupled system of partial differential equations (PDEs) subject to homogeneous Neumann boundary conditions. The system incorporates nutrient diffusion, fish chemotaxis toward nutrient-rich regions, and the filter-feeding role of mussels in nutrient recycling. The governing equations (1.1) include a parameter $\tau \in \{0, 1\}$, distinguishing between a parabolic-parabolic-parabolic system ($\tau=1$) and a parabolic-elliptic-parabolic system ($\tau=0$). The study aims to establish the global existence and boundedness of classical solutions, analyze conditions leading to finite-time blow-up (particularly in three dimensions), and develop a numerical scheme to validate these theoretical findings.

Methodology
The authors employ a combination of analytical techniques and numerical simulations:

1. Analytical Framework:

   • Local Existence: The existence of a unique local-in-time classical solution is established using semigroup methods and the Banach fixed-point theorem, relying on elliptic and parabolic regularity theory.
   • A Priori Estimates: To prove global existence, the authors derive $L^\infty$ and $L^p$ bounds. For the case $\tau=1$, a specialized test function $\phi(x)$ (constructed via Lemma 2.5) is used to handle the chemotactic term $-\chi\nabla \cdot (v\nabla u)$. This allows for the derivation of $L^p$ estimates for the fish population $v$ under specific constraints on the chemotactic sensitivity $\chi$ and nutrient bounds.
   • Extensibility Criteria: The proof of global existence relies on showing that the $L^\infty$ norm of the solution remains bounded, preventing blow-up. This is achieved by combining the $L^p$ estimates with the smoothing properties of the Neumann heat semigroup.
   • Blow-Up Analysis: For the three-dimensional case ($n=3$), the authors derive a lower bound for the blow-up time $T_{max}$. This involves constructing an auxiliary energy functional $\Psi_\tau$ (defined in 1.8) and deriving a differential inequality for its evolution. By integrating this inequality, an explicit lower bound for the blow-up time is obtained.

2. Numerical Approach:

   • Discretization: A finite element method (FEM) is developed using piecewise linear ($P1$) tetrahedral elements for spatial discretization and the backward Euler scheme for temporal integration.
   • Implementation: The fully coupled nonlinear system is solved using the PETSc SNES Newton line-search solver, with linearized systems handled by GMRES preconditioned with LU factorization. The simulations are performed using FEniCSx and visualized in MATLAB.
   • Validation: The numerical scheme is validated through mesh convergence studies and by comparing long-time dynamics against a spatially homogeneous ODE model solved via MATLAB's ode45.

Key Contributions and Results

• Global Existence and Boundedness (Theorem 1.1): The paper establishes the global existence of unique classical solutions for the system in dimensions $n \geq 3$. Specifically, if the initial data satisfies certain regularity conditions and the parameters satisfy the inequality:
  $$\max \left\{ \tau\|u_0\|_{L^\infty(\Omega)}, \frac{M}{\alpha} \right\} < \min \left\{ \frac{\pi}{\chi}\sqrt{\frac{2}{n}}, \frac{\delta_2}{\sigma_2} \right\}$$
  then the solution $(u, v, w)$ is uniformly bounded in time. This result holds for both $\tau=0$ and $\tau=1$.

• Blow-Up Time Estimates (Theorem 1.2): In the case where solutions blow up in finite time (specifically in a 3D convex domain), the authors derive a computable lower bound for the blow-up time $T_{max}$. This bound is expressed as an integral involving the initial value of the energy functional $\Psi_\tau(0)$ and system parameters. An explicit simplified lower bound is also provided:
  $$T_{max} \geq \frac{1}{2E_\tau \Psi_\tau(0)^2}$$

• Numerical Convergence: The proposed finite element scheme is shown to achieve second-order convergence in both $L^2(\Omega)$ and $L^\infty(\Omega)$ norms, consistent with the use of linear elements.

• Numerical Simulations:

  • Stability: Simulations with moderate parameters demonstrate that the system converges to a coexistence steady state, validating the long-time stability predicted by the ODE analysis.
  • Blow-Up Behavior: Numerical experiments in a 3D cubic domain with high initial concentrations and specific parameter choices illustrate the rapid, uncontrolled growth of the fish population density $v$, confirming the theoretical possibility of finite-time blow-up. The simulations show that while $v$ blows up, the nutrient $u$ and mussel $w$ densities remain bounded, consistent with the derived theoretical bounds.

Significance and Claims
The paper claims to provide a rigorous mathematical foundation for understanding the stability and potential instability of aquaculture ecosystems involving chemotactic fish and filter-feeding mussels. By extending previous ODE-based models to a PDE framework, the work highlights the critical role of spatial heterogeneity and directed movement (chemotaxis) in ecosystem dynamics.

The authors position their work as a bridge between theoretical analysis and computational modeling. They demonstrate that while global existence can be guaranteed under specific parameter regimes (preventing eutrophication-induced collapse), the system is susceptible to finite-time blow-up in higher dimensions if these conditions are violated. The development of a robust finite element scheme allows for the visualization of these complex dynamics, offering a tool to verify theoretical blow-up criteria and simulate the transition from stable coexistence to population explosion. The study does not propose new ecological management strategies but rather provides the analytical and numerical tools necessary to assess the stability of such systems under varying environmental conditions.