Sea-Cage Aquaculture, Sea Lice, and Declines of Wild Fish

Abstract: A sea cage, sometimes referred to as a net pen, is an enclosure designed to prevent farm fish from escaping and to protect them from large predators, while allowing a free flow of water through the cage to carry away waste. Farm fish thus share water with wild fish, which enables transmission of parasites, such as sea lice, from wild to farm and farm to wild fishes. Sea lice epidemics, together with recently documented population‐level declines of wild salmon in areas of sea‐cage farming, are a reminder that sea‐cage aquaculture is fundamentally different from terrestrial animal culture. The difference is that sea cages protect farm fish from the usual pathogen‐control mechanisms of nature, such as predators, but not from the pathogens themselves. A sea cage thus becomes an unintended pathogen factory. Basic physical theory explains why sea‐cage aquaculture causes sea lice on sympatric wild fish to increase and why increased lice burdens cause wild fish to decline, with extirpation as a real possibility. Theory is important to this issue because slow declines of wild fish can be difficult to detect amid large fluctuations from other causes. The important theoretical concepts are equilibrium, host‐density effect, reservoir‐host effect, and critical stocking level of farmed fish (stocking level at which lice proliferate on farm fish even if wild fish are not present to infect them). I explored these concepts and their implications without mathematics through examples from salmon farming. I also considered whether the lice‐control techniques used by sea‐cage farmers (medication and shortened grow‐out times) are capable of protecting wild fish. Elementary probability showed that (where W is the abundance of wild fish, W* is the prefarm abundance, F is the abundance of farm fish, and is the ratio of lice per farm fish to lice per wild fish). Declines of wild fish can be reduced by short growing cycles for farm fish, medicating farm fish, and keeping farm stocking levels low. Declines can be avoided only by ensuring that wild fish do not share water with farmed fish, either by locating sea cages very far from wild fish or through the use of closed‐containment aquaculture systems. These principles are likely to govern any aquaculture system where cage‐protected farm hosts and sympatric wild hosts have a common parasite with a direct life cycle.     

Characterization of bedload intermittency near the threshold of motion using a Lagrangian sediment transport model

At the smallest scales of sediment transport in rivers, the coherent structures of the turbulent boundary layer constitute the fundamental mechanisms of bedload transport, locally increasing the instantaneous hydrodynamic forces acting on sediment particles, and mobilizing them downstream. Near the critical threshold for initiating sediment motion, the interactions of the particles with these unsteady coherent structures and with other sediment grains, produce localized transport events with brief episodes of collective motion occurring due to the near-bed velocity fluctuations. Simulations of these flows pose a significant challenge for numerical models aimed at capturing the physical processes and complex non-linear interactions that generate highly intermittent and self-similar bedload transport fluxes. In this investigation we carry out direct numerical simulations of the flow in a rectangular flat-bed channel, at a Reynolds number equal to Re = 3632, coupled with the discrete element method to simulate the dynamics of spherical particles near the bed. We perform two-way coupled Lagrangian simulations of 48,510 sediment particles, with 4851 fixed particles to account for bed roughness. Our simulations consider a total of eight different values of the non-dimensional Shields parameter to study the evolution of transport statistics. From the trajectory and velocity of each sediment particle, we compute the changes in the probability distribution functions of velocities, bed activity, and jump lengths as the Shields number increases. For the lower shear stresses, the intermittency of the global bedload transport flux is described by computing the singularity or multifr actal spectrum of transport, which also characterizes the widespread range of transport event magnitudes. These findings can help to identify the mechanisms of sediment transport at the particle scale. The statistical analysis can also be used as an ingredient to develop larger, upscaled models for predicting mean transport rates, considering the variability of entrainment and deposition that characterizes the transport near the threshold of motion.     

Modeling and experiments of polydisperse particle clouds

A model for polydisperse particle clouds has been developed in this study. We extended the monodisperse particle cloud model of Lai et al. (Environ Fluid Mech 13(5):435–463, 2013) to the case of polydisperse particles. The particle cloud is first considered to be a thermal or buoyant vortex ring, with the thermal induced velocity field modeled by an expanding spherical Hill’s vortex. The buoyancy of the composite thermal is assumed to be the sum of buoyancy contributed by the all particles inside the thermal. Individual particles (of different particle properties) in the cloud are then tracked by the particle tracking equation using the computed induced velocity field. The turbulent dispersion effect is also accounted for by using a random walk model. Experiments of polydisperse particle clouds were carried out to validate the model. The agreement between model predictions and experiments was reasonable. We further validate our model by comparing it with the LES study of Wang et al. (J Hydraul Eng ASCE 141(7):06015006, 2014). The limitations of our model are then discussed with reference to the comparison. Overall, although some flow details are not captured by our model, the simplicity and generality of the model makes it useful in engineering applications.     

Kelvin–Helmholtz instability in the presence of variable viscosity for mudflow resuspension in estuaries

The temporal stability of a parallel shear flow of miscible fluid layers of different density and viscosity is investigated through a linear stability analysis and direct numerical simulations. The geometry and rheology of this Newtonian fluid mixing can be viewed as a simplified model of the behavior of mudflow at the bottom of estuaries for suspension studies. In this study, focus is on the stability and transition to turbulence of an initially laminar configuration. A parametric analysis is performed by varying the values of three control parameters, namely the viscosity ratio, the Richardson and Reynolds numbers, in the case of initially identical thickness of the velocity, density and viscosity profiles. The range of parameters has been chosen so as to mimic a wide variety of real configurations. This study shows that the Kelvin–Helmholtz instability is controlled by the local Reynolds and Richardson numbers of the inflection point. In addition, at moderate Reynolds number, viscosity stratification has a strong influence on the onset of instability, the latter being enhanced at high viscosity ratio, while at high Reynolds number, the influence is less pronounced. In all cases, we show that the thickness of the mixing layer (and thus resuspension) is increased by high viscosity stratification, in particular during the non-linear development of the instability and especially pairing processes. This study suggests that mud viscosity has to be taken into account for resuspension parameterizations because of its impact on the inflection point Reynolds number and the viscosity ratio, which are key parameters for shear instabilities.     

Testing Accuracy of Finite Element and Random Walk Schemes in Prediction of Pollutant Dispersion in Coastal Waters

Distribution of pollutants in coastal waters is usually represented by depth averaged twodimensional convection-dispersion equation. Under very specific conditions this equation can be solved analytically. Although such a solution is restricted to simplified situations it provides a very useful case for testing the performance of various numerical solution techniques currently available for the simulation of convective-dispersion of pollutants in natural water systems. In this paper the analytical solution of the convective dispersion equation is used as a benchmark against which the accuracy of other techniques are assessed. These assessments are based on quantitative comparisons between the results of the solution of two-dimensional convection-dispersion equation by the deterministic finite element and stochastic random walk methods. Both Eulerian and Lagrangian frameworks are employed to obtain the finite element solution of the convection-dispersion problem. It has been shown that the Lagrange–Galerkin finite element scheme yields the most accurate results for the case under study. However, computational costs of the Lagrange–Galerkin method can be relatively high and under certain conditions it may be reasonable to use a less accurate but cost effective random walk scheme to make water quality management decisions.     

An equation of motion for particles of finite Reynolds number and size

In order to simulate the motion of bubbles, drops, and particles, it is often important to consider finite Reynolds number effects on drag, lift, torque, and history force. Herein, an equation of motion is developed for spherical particles with a no-slip surface based on theoretical analysis, experimental data, and surface-resolved simulations. The equation of motion is then extended to account for finite particle size. This extension is critical for particles which will have a size significantly larger than the grid cell size, particularly important for bubbles, and low-density particles. The extension to finite particle size is accomplished through spatial-averaging (both volume-based and surface-based) of the continuous flow properties. This averaging is consistent with the Faxen limit for solid spheres at small Reynolds numbers and added mass and fluid stress forces at inviscid limits. The finite Re _p corrections are shown to have good agreements with experiments and resolved-surface simulations. The finite size corrections are generally fourth-order accurate and an order of magnitude more accurate than point-force expressions (which neglect quadratic and higher spatial gradients) for particles with size on the order of the gradient length-scales. However, further work is needed for more quantitative assessment of the truncation terms and the overall model robustness and accuracy in complex flows.     

Harvesting in an eight-species ecosystem

The theory for a general equilibrium ecosystem model that can include large number of interacting species is presented. Features include: (1) individual plants and animals are assumed to behave as if they are maximizing their net energy intake, (2) short- and long-run equilibriums are obtained, (3) species’ population adjustments depend on individual net energies. The theory is applied using simulations of an eight-species Alaskan marine ecosystem for which a “natural” equilibrium is calculated. Humans are introduced by adding a regulated open access fishery that harvests one of the species. Fishing impacts the fish population as well as the populations of other species, including Stellar sea lions, an endangered species. The sensitivity of fish and nonfish species populations to harvesting are calculated.     

Analytical expression for predicting the reduced settling velocity of small particles in turbulence

The influence of turbulence on the settling velocity of small particles remains an inconclusive research subject. Both enhanced and retarded particle settling compared to quiescent settling have been reported in previous literature, and several theories have been proposed. Among the mechanisms that account for reduced settling velocities in turbulence, the loitering effect is an important one that generally exists in various conditions. This study focuses on an analytical prediction of reduced particle settling velocities due to the loitering effect. By considering the velocity autocorrelation function as a step function analogous to the free path theory, and using a modified integral time scale to capture the essence of the loitering effect, an analytical expression is derived for predicting the reduction in particle settling velocities in turbulent flows. Calculation results of the expression are then examined by comparing with the results of a random walk model, direct numerical simulations that conditionally captured the reduction of particle settling velocities due to the loitering effect, and several representative experiments. Major possible influencing factors on the prediction of the analytical expression and applicable conditions of the expression are then further discussed. The proposed analytical expression is shown to be suitable for predicting the reduced settling velocities of small particles with relatively weak inertia in turbulent flows and could provide a reasonable explanation for reported cases in which small to moderate reductions in particle settling velocities were observed.

     

Part I : <Emphasis Type="Italic">a priori</Emphasis> study of erosion and deposition with large eddy simulation of turbulent flow over multiple 2D sandy Gaussian hills

The micro-scale prediction of sand trapping or take-off over hilly terrains is a crucial issue in semi-arid regions for soil depletion. In this context, large eddy simulations around one or several hills are performed in order to provide statistical parameters to characterize the flow at micro-scales and provide data for mesoscale modelling. We focus on the determination of recirculation zones since they play an important role in solid particle erosion or entrapment. A new wall modeling adapted from Huang et al. (J Turbul 17:1–24, 2016) for rough boundary layers is found to improve the prediction of the recirculation zone length downstream of an isolated hill and is used for all the numerical cases presented here. A geometrical parameterization of the recirculation zones is proposed. When the recirculation region is assumed to have an ellipsoidal shape, the total surface of the recirculation can be obtained from this new parameterization and easily extrapolated to more general dune configurations. Numerical results are compared with experiments performed in our laboratory (Simoëns et al. in Procedia IUTAM 17:110–118, 2015) and good agreement is achieved. We explore general aerodynamic cases deduced from the urban canopy scheme of Oke (Energy Build 11:103–113, 1988). In this scheme the momentum and mass exchange between the upper layer and the space between hills is sorted according to the streamwise hill spacing within three basic cases of skimming, wake or isolated flow. The study of the recirculation zones, the mean velocity and Reynolds stress profiles around an isolated or two consecutive hills with different distances shows that the double hill configuration with 3H separation behaves as much as a whole to the upcoming flow. The vortex formed between the crests does not strongly affect the overall evolution of the outer flow. By an a priori prediction of the preferential zones of erosion and accumulation of fictive particles, it is shown that isolated dunes present more deposition and less erosion than two-hill configurations. The results presented in this study will be discussed in the presence of Lagrangian transport of sand particles above 2D Gaussian hills in future work.     

A general random walk model for the leptokurtic distribution of organism movement: Theory and application

The movement of organisms is usually leptokurtic in which some individuals move long distances while the majority remains at or near the area they are released. There has been extensive research into the origin of such leptokurtic movement, but one important aspect that has been overlooked is that the foraging behaviour of most organisms is not Brownian as assumed in most existing models. In this paper we show that such non-Brownian foraging indeed gives rise to leptokurtic distribution. We first present a general random walk model to describe the organism movement by breaking the foraging of each individual into events of active movement and inactive stationary period; its foraging behaviour is therefore fully characterized by a joint probability of how far the individual can move in each active movement and the duration it remains stationary between two consecutive movements. The spatio-temporal distribution of the organism can be described by a generalized partial differential equation, and the leptokurtic distribution is a special case when the stationary period is not exponentially distributed. Empirical observations of some organisms living in different habitats indicated that their rest time shows a power-law distribution, and we speculate that this is general for other organisms. This leads to a fractional diffusion equation with three parameters to characterize the distributions of stationary period and movement distance. A method to estimate the parameters from empirical data is given, and we apply the model to simulate the movement of two organisms living in different habitats: a stream fish (Cyprinidae: Nocomis leptocephalus) in water, and a root-feeding weevil, Sitona lepidus in the soil. Comparison of the simulations with the measured data shows close agreement. This has an important implication in ecology that the leptokurtic distribution observed at population level does not necessarily mean population heterogeneity as most existing models suggested, in which the population consists of different phenotypes; instead, a homogeneous population moving in homogeneous habitat can also lead to leptokurtic distribution.     

Zooplankton encounters in patchy particle distributions

Zooplankton encounter rates are dependent not only on both sensory and swimming performances of the organisms, but also on the distribution pattern of food particles. Increasing evidences indicate that, in natural conditions, phytoplankton is often aggregated in thin layers. In the present contribution we investigate the concomitant effects of motion complexity and habitat fragmentation on the number of encounters realised by virtual continuously moving copepods adopting different motion strategies. Our simulated organisms move in an environment characterised by the presence of thin patches of phytoplankton, and their swim follows five motion rules (pure random walk, correlated random walk with three different time scales, self-avoiding random walk), each characterised by a typical value of the three-dimensional fractal dimension D_3D. Compared to a uniform distribution, for a given motion rule the clustering of prey particles increases the variance of encounters, while no remarkable effect is reported in the average number of particles intercepted. These results broaden our understanding of the behavioural efficiency in freely swimming zooplankters and improve our knowledge of the functioning of aquatic systems.     

The evolution of large and small-scale structures in Kelvin–Helmholtz instabilities

To better understand the dynamics of Kelvin–Helmholtz instabilities in environmental flows, their evolution is investigated using direct numerical simulations (DNS). Two-dimensional DNS is used to examine the large-scale and small-scale structures of the instability at high Reynolds and Prandtl numbers that represent real environmental flows. The semi-analytical model of Corcos and Sherman (J Fluid Mech 73:241–264, 1976) is used to explain the physics of these simulations prior to saturation of the KH billow, and also provide a computationally efficient prediction of the vortex dynamics of the instability. The DNS results show that the large-scale structure of the billow does not depend on the Reynolds number for sufficiently high Reynolds numbers. The billow structure reveals a less straightforward dependence on the Prandtl number. Predictions of the model of Corcos and Sherman (J Fluid Mech 73:241–264, 1976) improve as Reynolds number and Prandtl number increase. The small-scale structure of the vorticity and density fields vary with both Reynolds and Prandtl numbers. Three-dimensional DNS of KH flows and their transition to turbulence are used to study small length scales. Based on the thickness of the braid, a simple method is introduced to estimate the Batchelor scale, which can be used as a guide for the resolution required for the direct numerical simulation of two and three-dimensional Kelvin–Helmholtz flow fields.     

Vertical fiscal externalities and the environment

We show that imposition of a state-level environmental tax in a federation crowds out pre-existing federal taxes. We explain how this vertical fiscal externality can lead unilateral state-level environmental policy to generate a welfare gain in the implementing state, at the expense of other states, even absent any environmental benefits. Using a computable general equilibrium model of the Canadian federation, we show that vertical fiscal externalities can be the major determinant of the welfare change following environmental policy implementation by a state government. Our numerical simulations indicate that – as a consequence of vertical fiscal externalities – state governments can reduce greenhouse gas emissions by over 20 percent without any net cost to themselves.     

Mixing at a sediment concentration interface in turbulent open channel flow

In this work we address the role of turbulence on mixing of clear layer of fluid with sediment-laden layer of fluid at a sediment concentration interface. This process can be conceived as the entrainment of sediment-free fluid into the sediment-laden layer, or alternatively, as the transport of sediment into the top sediment-free flow. This process is governed by four parameters—Reynolds number of the flow \(Re_\tau\), non-dimensional settling velocity of the sediment (proxy for sediment size) \(\tilde{V}\), Richardson number \(Ri_\tau\) and Schmidt number Sc. For this work we have performed direct numerical simulations for fixed Reynolds and Schmidt numbers while varying the values of Richardson number and particle settling velocity. In the simple model considered here, the flow’s momentum and turbulence pre-exists over the entire layer of fluid, while the sediment is initially confined to a layer close to the bed. Mixing of sediment-free fluid with the sediment-laden layer is associated primarily with upward transport of sediment and buoyancy. There is no simultaneous upward transport of fluid momentum and turbulence into the sediment-free fluid layer, which is already in motion and turbulent. The analysis performed shows that the ability of the flow to transport a given sediment size decreases with the distance from the bottom, and thus only fine enough sediment particles are transported across the sediment concentration interface. For these cases, the concentration profiles evolve to a final steady state in good agreement with the well-known Rouse profile. The approach towards the Rouse profile happens through a transient self-similar state. This behavior of the flow is not seen for larger particles. Detailed analysis of the three dimensional structure of the sediment concentration interface shows the mechanisms by which sediment particles are lifted up by tongues of sediment-laden fluid with positive correlation between vertical velocity and sediment concentration. Finally, the mixing ability of the flow is addressed by monitoring the time evolution of the center of mass of the sediment-laden layer and the vertical location of the sediment-free/sediment-laden interface.     

Wind- and tide-induced currents in the Stagnone lagoon (Sicily)

The hydrodynamic circulation is analyzed in the coastal lagoon of Stagnone di Marsala, a natural reserve located in the north-western part of Sicily, using both experimental measurements and numerical simulations. Field measurements of velocities and water levels, carried out using an ultrasound sensor (3D), are used to validate the numerical model. A 3D finite-volume model is used to solve the Reynolds-averaged momentum and mass balance differential equations on a curvilinear structured grid, employing the k–e{\varepsilon} turbulence model for the Reynolds stresses. The numerical analysis allows to identify the relative contribution of the forces affecting the hydrodynamic circulation inside the lagoon. In the simulations only wind and tide forces are considered, neglecting the effects of water density changes. Two different conditions are considered. In the first both the wind stress over the free-surface and the tidal motion are imposed. In the second the wind action is neglected, to separately analyze the tide-induced circulation. The comparison between the two test cases highlights the fundamental role of the wind on the hydrodynamics of the Stagnone lagoon, producing a strong vertical recirculation pattern that is not observed when the flow is driven by tides only.     

Incorporating behaviour into simple models of dispersal using the biological control agent Dicyphus hesperus

We explored the utility of incorporating easily measured, biologically realistic movement rules into simple models of dispersal. We depart from traditional random walk models by designing an individual-based simulation model where we decompose animal movement into three separate processes: emigration, between-patch movement, and immigration behaviour. These processes were quantified using experiments on the omnivorous insect Dicyphus hesperus moving through a tomato greenhouse. We compare the predictions of the individual-based model, along with a series of biased random walk models, against an independent experimental release of D. hesperus. We find that in this system, the short-term dispersal of these insects is described well by our individual-based model, but can also be described by a 2D grid-based biased random walk model when mortality is accounted for.     

A Computational Fluid Dynamics Approach for Urban Area Transport and Dispersion Modeling

A simulation tool has been developed to model the wind fields, turbulence fields, and the dispersion of Chemical, Biological, Radiological and Nuclear (CBRN) substances in urban areas on the building to city blocks scale. A Computational Fluid Dynamics (CFD) approach has been taken that naturally accounts for critical flow and dispersion processes in urban areas, such as channeling, lofting, vertical mixing and turbulence, by solving the steady-state, Reynolds-Averaged Navier–Stokes (RANS) equations. Rapid generation of high quality cityscape volume meshes is attained by a unique voxel-based model generator that directly interfaces with common Geographic Information Systems (GIS) file formats. The flow and turbulence fields are obtained by solving the steady-state RANS equations using a collocated, pressure-based approach formulated for unstructured and polyhedral mesh elements. Turbulence modeling is based upon the Renormalization Group variant of the k–ε model (k–ε RNG). Neutrally buoyant simulations are made by prescribing velocity boundary condition profiles found by a power–law relationship, while turbulence quantities boundary conditions are defined by a prescribed mixing length in conjunction with the assumption of turbulence equilibrium. Dispersion fields are computed by solving an unsteady transport equation of a dilute gas, formulated in a Eulerian framework, using the velocity and turbulence fields found from the steady-state RANS solution. In this paper the model is explained and detailed comparisons of predicted to experimentally obtained velocity, turbulence and dispersion fields are made to neutrally stable wind tunnel and hydraulic flume experiments.     

Diffusion about the mean drift location in a biased random walk

Random walks are used to model movement in a wide variety of contexts: from the movement of cells undergoing chemotaxis to the migration of animals. In a two-dimensional biased random walk, the diffusion about the mean drift position is entirely dependent on the moments of the angular distribution used to determine the movement direction at each step. Here we consider biased random walks using several different angular distributions and derive expressions for the diffusion coefficients in each direction based on either a fixed or variable movement speed, and we use these to generate a probability density function for the long-time spatial distribution. We demonstrate how diffusion is typically anisotropic around the mean drift position and illustrate these theoretical results using computer simulations. We relate these results to earlier studies of swimming microorganisms and explain how the results can be generalized to other types of animal movement.     

Froude scaling limitations in modeling of turbidity currents

The scaling problem associated with the modeling of turbidity currents has been recognized but is yet to be explored systematically. This paper presents an analysis of the dimensionless governing equations of turbidity currents to investigate the scale effect. Three types of flow conditions are considered: (i) conservative density current; (ii) purely depositional turbidity current; and (iii) mixed erosional/depositional turbidity current. Two controlling dimensionless numbers, the Froude number and the Reynolds number, appear in the non-dimensional governing equations. When densimetric Froude similarity is satisfied, the analysis shows that the results would be scale-invariant for conservative density current under the rough turbulent condition. In the case of purely depositional flows, truly scale-invariant results cannot be obtained, as the Reynolds-mediated scale effects appear in the bottom boundary conditions of the flow velocity and sediment fall velocity. However, the scale effect would be relatively modest. The Reynolds effect becomes more significant for erosional or mixed erosional/depositional turbidity currents as Reynolds-mediated scale effects also appear in the sediment entrainment relation. Numerical simulations have been conducted at three different scales by considering densimetric Froude scaling alone as well as combined densimetric Froude and Reynolds similarity. Simulation results confirm that although the scaling of densimetric Froude number alone can produce scale-invariable results for conservative density currents, variations occur in the case of turbidity currents. The results become scale invariant when densimetric Froude and Reynolds similarities are satisfied simultaneously.     

Spatial memory, habitat auto-facilitation and the emergence of fractal home range patterns

Animals interact with their habitat in a manner which involves both negative and positive feedback mechanisms. We apply a specific modeling approach, “multi-scaled random walk”, for the scenario where a spatially explicit positive feedback process emerges from a combination of a spatial memory-dependent tendency to return to familiar patches and a consequently objective or subjective improvement of the quality of these patches (habitat auto-facilitation). In addition to the potential for local resource improvement from physically altering a patch, primarily known from the ecology of grazing ungulates, auto-facilitation from site fidelity may also embed more subtle subjective, individual-specific advantages from patch familiarity. Under the condition of resource superabundance, fitness gain from intra-home range patch fidelity creates a self-reinforcing use of the preferred patches on expense of a broader foraging in a priori equally favorable patches. Through this process, our simulations show that a spatially fractal dispersion of accumulated locations of the individual will emerge under the given model assumptions. Based on a conjecture that intra-home range patch fidelity depends on spatial memory we apply the multi-scaled random walk model to construct a spatially explicit habitat suitability parameter H_ij, which quantifies the dispersion of the generally most constraining resource from the individual's perspective. An intra-home range set of observed H-scores, H_obs, can then be estimated from a simple 2-scale calculation that is derived from the local dispersion of fixes. We show how the spatially explicit habitat utilization index H_obs not necessarily correlates positively with the local density fluctuations of fixes. The H-index solves some well-known problems from using the pattern of local densities of telemetry fixes - the classic utilization distribution - as a proxy variable for relative intra-home range habitat quality and resource selection. A pilot study on a set of telemetry fixes collected from a herd of free-ranging domestic sheep with overlapping summer home ranges illustrates how the H-index may be estimated and interpreted as a first-level approach towards a more extensive analysis of intra-home range habitat resource availability and patch preferences. Spatial memory in combination with site fidelity requires a modeling framework that explicitly describes the property of positive feedback mechanism under auto-facilitation in a spatio-temporally explicit manner.