An Ordinary Differential Equation (ODE) that governs the longitudinal movement of a fish school is presented. The 1-D domain D = (0, L) with the length L > 0 is considered as a river reach bounded by physical barriers placed at the boundary ∂D that contains the upstream- and downstream-ends x = 0, L. The flow velocity V:D¯→ℝ is positive and Lipschitz continuous in D¯. The fish school is considered as a synchronized group, which is seen as a moving point in D. The position of the school at the time t ≥ 0 is denoted as X_t: [0, T]→ D¯ where the time that X_t firstly hits ∂D is denoted as T. The swimming velocity u_t: [0, T] → U = [-u_max, u_max] of the fish school is considered as a control variable. Here, u_max is the maximum sustained swimming speed, which is assumed to be sufficiently large such that umax>maxx∈D¯V(x). The ODE that governs the longitudinal movement of the fish school until it first hits the boundary ∂D is set as:

The performance indexed to be maximized by the fish school, the decision-maker, is presented. It is assumed that all the individuals in the school share a common performance index containing the three terms: benefit of finding a high-quality habitat, hydrodynamic cost, and non-hydrodynamic cost. The population density in the fish school N_t : [0, T] → [0, +∞], referred to as the size in what follows, is a function of the time t. The stopping time of migration is denoted as τ. In what follows, α∈C1(D¯ ) represents the profit by terminating the migration, which is assumed to be non-negative. The performance index J is set as:

with the discount rate δ ≥ 0. Conceptually, δ represents impatience of the fish school: larger δ results in earlier stopping of migration. It would be natural to consider that the existence of a performance index to be maximized during a fish migration process is a result of evolutionary processes. Also, δ can be larger for the migration subject to higher predation pressure: namely higher mortality.

In eq. (2), the function f ≥ 0 with f(0) = 0 is the hydraulic cost of solo swimming per unit time as a smooth, non-negative, convex, and even function. This f can be identified from laboratory or field experiments on the swimming behavior of individual fishes along a current as demonstrated later. The coefficient a > 0 represents the discount of the hydrodynamic cost by forming a school, which is motivated by the theoretical and experimental results that increasing the size N effectively reduces the hydrodynamic cost per individual [23]. The probabilistic model theoretically predicts [24]:

On the other hand, the coefficient b > 0 conceptually represents the non-hydraulic cost by forming a school per unit time. Schooling would have many aspects of benefits, such as improvement of navigational performance, hearing perception, and foraging efficiency. While at the same time, schooling would negatively affect passage efficiency, information transmission among individuals, and competitions among the individuals [25]. Being different from the sophisticated multi-agent models [26], a simple strategy is employed in this paper. This is because these effects are possibly interacting and difficult to mathematically describe in detail. Therefore, the lumped approach [27] is employed:

with k > 0 where its sophisticated parameterization remains as a future topic. This term is necessary so that a fish school is created. For example, removing this term from the performance index J leads to an unreasonable result that forming a fish school with the infinitely large fish school is optimal.

The value function Φ:D¯→ℝ is the maximized performance index defined as:

From a biological viewpoint, the value function is an index to be maximized by the fish school through choosing a migration strategy. The optimizers (u, N, τ) to give the supremum of eq. (5) are expressed as (u_opt, N_opt, τ_opt), which are functions of x. Finding the optimizers is the goal of the present model.

The dynamic programming principle (Chapter III.4 of [18]) leads to the VI as a governing differential equation of the value function Φ:

where H:D¯×ℝ×ℝ→ℝ is the function called Hamiltonian, and the function F:D¯×ℝ→ℝ is given as:

The derivative of Φ with respect to x is denoted as Φ'. Hereafter, the function f is assumed to lead to the inequality:

where c₁, c₂ > 0 are some constants depending on neither x nor p. Also, the condition k > m is assumed to make convex be (fa⁻¹ + b) (u, N). The latter assumption leads to a concave F(x, p) for p ∈ ℝ , and the equation F(x, p) has the two solutions p = p⁻ (x), 0 with p⁻(x) < 0 for x e ∈ D¯. At least two examples of f that comply with the conditions mentioned above are found, both of which are used in this paper: f1(u) = |u|ⁿ⁺¹ with kn > m (condition of forming a school) and f2(u)=(1−1−|u|/umax)2 [27]. eq. (8) is a key for the unique solvability of eq. (6).

Viscosity solutions are appropriate candidates of weak (non-classical) solutions to degenerate elliptic differential equations like the eq. (6). This VI does not have classical, continuously differentiable solutions in general but has viscosity solutions that are not differentiable and possibly discontinuous. Handling the eq. (6) requires the concept of viscosity solutions even under simplified conditions as shown in the next sub-section. For the sake of brevity of analysis, the boundary operator B : ∂D×ℝ → ℝ is introduced as:

moreover, the operators G∗,G∗:D¯×ℝ×ℝ→ℝ as:

and:

for the sake of brevity of descriptions. The upper- and lower semi-continuous envelops of Φ : D¯ → ℝ are denoted as Φ^* and Φ*, respectively [28]. They are defined as:

and:

respectively, where Β_ε(x) with ε > 0 represents the closed interval [x − ε, x + ε].

Following the literature for viscosity solutions [27], the definition of viscosity solutions to the eq. (6) is stated as follows.

An exact viscosity solution to the eq. (6) is derived for specific functional forms of the coefficients, which can give biological and ecological implications of the proposed mathematical model. Assume the uniform flow condition V = const < u_max and:

with the constants A, B, C > 0. This α is motivated by the fact that the quality of habitat measured by the amount of food, such as diatoms, would sharply increase along a river [29]. A formal exact solution to the eq. (6) on D¯ is then found as Figures 1a and 1b:

and:

Figure 1.

Exact (blue) and numerical solutions with the cell size l = 0.01 L (green), and α (black) for: continuous (A = 1) (a) and discontinuous case (A = 2) with f = f₁ and n = 2 (b); the other parameters are V = 1, u_max = 5, δ = 0, m = 1/3, k = 0.5, d = 0.5, B = 10, C = 20, λ = 10⁶ (the numerical solutions are computed with 100 uniform cells)

In eq. (16), p⁻ is a negative constant determined as a solution of F = 0. The two graphs y = γ(x) and y = α(x) intersect at x = y₁ such that y₀ < y₁ ≤ L if:

Therefore, when eq. (18) is satisfied, the found Φ in eq. (16) is rewritten for x > y₀ as:

Similarly, when eq. (18) is not satisfied, Φ in eq. (16) is rewritten for x > y₀ as:

The optimal controls can be determined in the sub-domain of D¯ where α(x) < Φ(x). Assume that the minimum of the second term of eq. (7) is achieved by the internal solutions u = u_opt ∈ (−u_max, u_max ) and N = N_opt > 0. Then, they are positive constants and satisfy:

and:

To prove that the exact solution in eq. (16) is a viscosity solution is not a trivial issue because this Φ is not smooth at x = y₁. Since the exact solution satisfies the eq. (6) where the solution is continuously differentiable, it has to be examined with the condition of viscosity solution at which it is not differentiable or discontinuous. After elementary calculations, it turns out that eq. (16) is a viscosity solution to the eq. (6) by examining it against the condition of viscosity super- and sub-solutions on D¯.

Practical implications of the exact solution of eq. (16) are presented focusing on its dependence on the parameters and coefficients. The condition ω ≥ 1 implies that the cost of migration is sufficiently small and/or the benefit of finding a high-quality habitat is sufficiently high. The fish school approaches the upstream-end x = 0 when A is sufficiently high. The school does not migrate upstream (τ_opt = 0) for all x ∈ D¯ when ω < 1. Therefore, the exact solution in eq. (16) implies that the fish school would not migrate upstream along the river if its environment is severely degraded. In particular, y0 is increasing with respect to A, showing that degradation of the river environment results in more downstream termination of the fish migration. Therefore, the shape and the magnitude of the habitat quality α are important for qualitative understanding of the fish migration. This point is investigated in the fourth section as well through a model application.

For the cost function f₁ presented in the previous section, assuming the condition kn > m yields that u_opt and N_opt are increasing with respect to the flow speed V, which is consistent with the experimental results [30]. Similar statement holds for f₂. The second of eq. (21) states that the swimming speed is increasing with respect to the school size, which is also consistent with the experimental results [30].

In general, there is no guarantee that the found exact solution is the unique viscosity solution to the eq. (6). This is not only a mathematical issue but also a practical issue since the viscosity solution should be the plausible solution to the eq. (6). The issue of unique solvability is discussed below for the cases δ > 0 and δ = 0 separately. In this section, α of eq. (15) is employed.

For δ > 0, the eq. (6) is a unique continuous viscosity solution by the comparison theorem (Theorem 3.17 of [28]) since the eq. (6) satisfies the continuity and ellipticity conditions {eqs. (3.1), (3.24), and (3.35) of [28]}. The uniqueness of discontinuous viscosity solutions for δ > 0 are unclear at this stage, but it is expected that this subject can be dealt with based on the constructive argument presented below.

For δ = 0, the above comparison argument cannot be used since it requires δ > 0. This is the major difficulty for dealing with the case with δ = 0. The remaining part of this sub-section presents the idea and plan of a constructive proof of unique solvability of the eq. (6) for δ = 0 in the viscosity sense for the exactly-solvable problem in the previous sub-section. The constructive proof follows a step-by-step manner presented in Propositions 1 through 4. The proof is a bit long, but is presented in this paper for the sake of its self-containedness.

For the sake of brevity of descriptions, the sets D₁ and D₂ are defined as:

It is assumed that viscosity solutions do not have dense non-differentiable points in D. This assumption has a technical aspect, but viscosity solutions having dense non-differentiable points in D can be very irregular and not reasonable from biological viewpoints.

Solutions to the eq. (6) may not be explicitly found in applications, which motivates us to develop a numerical scheme that can approximate its solutions. For this purpose, a finite difference scheme based on exact viscosity solutions to local two-point boundary value problems is presented for solving the eq. (6) under general conditions. This kind of discretization technique utilizing local exact solutions is called fitting technique and has been applied to elliptic and parabolic problems [32]. The domain D is discretized into I ≥ 2 cells C_i = (x_i−1, x_i) (1 ≤ i ≤ I) with 0 = x₀ < x₁ < ... < x₁ = L. Φ is approximated at each x_i. The approximation of Φ at x_i is denoted as Φ_i. The length of C_i is denoted as li = x_i − x_i−1. The discretization procedure below focuses on a vertex x_i (1 ≤ i ≤ I − 1) without the loss of generality.

The finite difference scheme for a linear problem is firstly explained and is extended so that the present VI is handled. Consider the linear advection-decay equation:

with a piece-wise continuous g. W (≠ 0) and g are approximated in each cell and those in C_i are denoted with the subscript i as W_i and g_i, respectively. A finite difference approximation of eq. (24) is proposed as:

where F_i,j is expressed with Φ_i+j−1, Φ_i+j+1, W_i+j, and g_i+j. Consider the local linear two-point boundary value problem in C_i to find Ψ_i : Ψi:C¯ i→ℝ :

Its unique viscosity solution in C_i is found as:

The present scheme evaluates F_i,j with Ψ_i as:

The boundary conditions are directly specified at x₀ = 0 and x₁ = L. Assembling eq. (25) for 1 ≤ i ≤ I −1 with the boundary conditions leads to a system of linear equation with a coefficient matrix whose inverse is positive-definite, showing that the present scheme is unconditionally stable. The scheme reduces to the classical first-order upwind finite difference scheme as R → +0.

For solving the eq. (6), a penalty method [33] is employed to reduce a VI to a differential equation that is easier to numerically handle. The penalized counterpart of the eq. (6) is given as:

with the penalty parameter λ > 0, so that it is expressed like eq. (24). The original eq. (6) is formally recoveredwhen λ → +∞, which is implemented numerically in the penalty method. In the present scheme, the first and second terms of eq. (29) are advection-decay terms having solution-dependent coefficients. These coefficients are determined by optimizing the “min” term of eq. (7) in each cell. The third term is discretized at each vertex as λ max {αi − Φ_i, 0}. A standard fixed-point iteration is used for solving eq. (29). The iteration is terminated when the difference of old and updated numerical solutions at each vertex becomes smaller than 10⁻⁸, which is a sufficiently small value.

The present finite difference scheme is degenerate elliptic (Definition 2), Lipschitz continuous (Definition 3), and proper (Definition 6) in the sense of Oberman [34]. Therefore, for the linear problem of eq. (24), the scheme admits a unique numerical solution (Theorem 8 [34]). In addition, the scheme is monotone and non-expansive in the l^∞-norm since it is degenerate elliptic (Theorem 3 [34]). The same statements hold for eq. (29).

Convergence of the scheme for our VI is numerically examined with both continuous and discontinuous viscosity solutions. The domain is uniformly discretized into I cells. For the case δ = 0, the non-oscillatory numerical solutions presented in Figures 1a and 1b demonstrate its accuracy and stability. The numerical solutions accurately capture the non-differentiable and discontinuous points of the solutions. Tables 1 and 2 summarize the computed l^∞ errors for a series of λ and I⁻¹. In Table 1, “Osc” means non-convergence of numerical solutions. This also applies to the other tables in this paper. The computational results suggest the efficient scaling relationship λ = O(I⁻¹), which may not be a sharp relationship but useful for accurate numerical computation. Rigorous proof for this scaling result should be analyzed in future research.

For the case δ > 0, since our VI does not admit exact solutions, the scheme is examined against the following test problem with constant coefficients:

which can also be cast into ψ=1−VRL−1AR a penalized form as in eq. (29). The viscosity solution to eq. (30) is:

which is discontinuous at x = L when Ψ < 1, where:

In the numerical computation, the parameter values employed are V = 1, R = 2, L = 1, and A = 1 for the smooth, classical solution case and A = 10 for the non-smooth, viscosity solution case. Tables 3 and 4 summarize the computed l^∞ errors for a series of λ and I⁻¹. The computational results presented in this sub-section indicate the satisfactory performance of the present finite difference scheme. The error for the smooth, classical solution case is negligible. This is owing to the specialized discretization employed in the present finite difference scheme that the solutions with the exponential profiles like that in eq. (32) can be exactly reproduced.

The present mathematical model and the verified finite difference scheme are applied to the evaluation of upstream fish migration in Hii River, San-in area, Japan (Figure 2, top left). This river provides indispensable water resources for human activities in this area. The total length of the mainstream and the catchment area of Hii River are 153 km and 2,070 km², respectively. The river flows into the brackish Lake Shinji, which has a downstream brackish lake named Lake Nakaumi. Hydrology of Hi River are reviewed in Somura et al. [35]. Both brackish lakes are key habitats for migrating aquatic species that grow and spawn in Hii River. Hii River Fisheries Cooperatives (HRFC) authorizes inland fishery of the middle to upper reaches of the river.

Figure 2.

Map of the Hii River (top left), photos of Obara Dam (top right) and Hinobori Weir (bottom) (photo credits: H. Yoshioka)

P. altivelis (Ayu) (Figure 3, top left) is one of the common diadromous fishes in Japan and is a major inland fishery resource in the country. The natural P. altivelis has an annual life cycle [36]. During autumn, adults spawn eggs in the downstream reaches of the living river and die soon afterwards. Hatched larvae descend to coastal areas of the downstreamwater body of the river, typically a sea or a brackish lake. The larvae grow up to juveniles with feeding on zoo planktons till the next spring. The grown fishes ascend the river toward its midstream and upstream reaches where diatoms (Figure 3, top right), which are staple foods of P. altivelis, are available on the riverbed. They feed on the algae to mature till the coming autumn when they descend the river.

Figure 3.

Photos of P. altivelis (top left), diatom (top right), and Cladophora glomerata Kützing (bottom) just downstream of Obara Dam (photo credits: H. Yoshioka)

Recently, fish catches of P. altivelis in Japan have been dramatically decreasing due to severe degradations of the river environment and ecosystems. Installing physical barriers, such as huge dams and weirs, into river cross-sections prevents fishes from upstream migration and significantly affects downstream water environment. The population of the fish in some rivers are maintained through intensive artificial hatching of farmed juveniles during the spring season in each year by local fishery cooperatives. The artificially hatched P. altivelis also migrate toward upstream after the hatching and grow till the coming autumn, but has been thought to be unsuccessful in reproduction due to genetic reasons. This would be the main reason why the fish have to be artificially hatched in a river in each spring to maintain their population.

A major issue in releasing the fish is when and where to release farmed P. altivelis along the river and its tributaries so that the released fishes grow well and the fish catch in the river increases. In Hii River, local fishery cooperatives and residents say that the harmful attached algae Cladophora glomerata Kützing (Figure 3, bottom) are significantly growing from the just downstream point toward several km downstream reach of the huge multi-purpose dam named Obara Dam (Figure 2, top right) after its construction in 2011. This dam does not have facilities like fishways that fishes can pass through, serving as a physical barrier. According to residents and the officers of HRFC, the areas at and around which the dam was created, and the harmful algae are growing involved good fishing ground of P. altivelis. They also report that the P. altivelis after the construction of the dam are significantly smaller than those before the construction. In addition, the experimental research indicated that P. altivelis could not digest the harmful algae [37]. The Hinobori Weir with the height of 11 m (Figure 2, bottom) for erosion control is the second largest physical barrier in the mainstream. The weir has a vertical-slot fishway longer than 100 m, but its passage efficiency is not clear. In each spring, many P. altivelis are artificially released into the river reach between the Hinobori Weir and Obara Dam. Currently, HRFC is faced with the problem when and where to artificially hatch farmed P. altivelis given the potential shift of the attached algae. This problem is partly approached in this paper.

The cost function f for P. altivelis is firstly presented, which is used in the simulation below. Yoshioka et al. [27] collected the observed datasets of swimming behavior of individual P. altivelis against water currents, and found the concave formula:

where V_flow ∈ [0, u_max] is the flow speed of the current and u_swim : [0, u_max] → [0, u_max] is the swimming velocity of the individual against the current with the speed V_flow. Correlation between the observed and modelled u_swim is 0.80 [27], showing its reasonable accuracy. The cost function f : [0, u_max] → [0, +∞] is then determined as a solution to the initial value problem of the ODE [27]:

which turns out to give f = c₃f₂(u) with a constant of integration c₃ that can be set as c₃ = 1 without loss of generality. f₂ is twice differentiable and convex in (0, u_max), and its derivative diverges at u = u_max. A limitation of the above u_swim and f is that they cannot simulate upstream fish migration with the flow speed greater than u_max. It is conjectured that, above the flow speed u_max, the fish is not able to maintain prolonged swimming speed and must perform burst swimming such that the present mathematical model cannot be directly used. What is important here is that the hydraulic cost function f is derived from the swimming behavior of the fish.

The river reach from just downstream of Obara Dam (x = 0 km) to just upstream of Hinobori Weir (x = L = 13 km) is set as the domain D, which is discretized into 218 cells with 219 nodes. A steady flow field in the Hii River system involving the river reach is computed with the 1-D shallow water solver [37] specifying the Manning’s friction coefficient of 0.05 s/m^1/3. This numerical solver has already been verified against the benchmark, experimental, and real cases. The flow discharge in the reach is 12.5 to 15.0 m³/s, which does not significantly deviate from the observed discharge during April to May. The flow velocity in the reach is utilized as the coefficient V for the present model (Figure 4). The profiles of 0 ≤ α ≤ αmax in D with the maximum value αmax > 0 are specified on the basis of the summation of two hyperbolic-tangential functions in eq. (15) whose coefficients were determined from the interviews from the officers and members of HRFC (Figure 4), α ≈ α_max at x∈ D¯ represents that diatoms are reported to be abundant at x, α ≈ 0 represents that diatoms are reported to be sparse at x, and otherwise the diatoms moderately exist. This more accurate identification of α will be carried out in our future research. The other model parameters are set as δ = 10⁻⁶ 1/s, m = 1/3, k = 2.0, d = 10⁻³ m²/s², and λ = 10⁶ 1/s.

Figure 4.

Computed flow velocity V and the modelled coefficient α along the reach

Figure 5 shows the computed area in the river where the fish school terminates its migration, namely the area where Φ = α in D¯ for a range of α_max. The results show that the area expands as α_max decreases. For moderately large α_max, there is a white area around 0.6L < x < 0.9L where the fish school ascends. The results show that hatched fishes in the downstream grey area around 0.9L < x < L would not ascend the river reach. In the other white area around 0 < x < 0.5L, the fish school swims toward the upstream-end of the potential habitat. In addition, a decrease of α_max, namely degradation of the overall quality of potential habitat of the fish as good fishing grounds, shows that the downstream-end of the potential habitat move toward upstream. Increasing the value of δ leads to the wider black area than that in Figure 5 as shown in Figure 6 for δ = 10⁻⁴ 1/s, and δ > 10⁻³ 1/s results in no upstream migration with Φ = α on D¯ except for near x = 0. Management of river environment is implied to be crucial for fish migration and inland fisheries.

Figure 5.

The computed area where the fish school terminates its migration along the river (lateral axis) with respect to α_max (vertical axis), the area is coloured blue, the parameter value is δ = 10⁻⁵ l/s in this figure

Figure 6.

The computed area where the fish school terminates its migration along the river (lateral axis) with respect to α_max (vertical axis), the area is coloured blue, the parameter value is δ = 10⁻⁴ 1/s in this figure

The computational results suggest the following management policy for the studied river reach. Especially for the upstream part of the river reach (0 < x < 0.5L), the habitat quality should be improved so that this part serves as a habitat for the fish. For example, cleaning up the riverbed to exterminate the harmful algae can be an effective way for that purpose. The results also suggest that the fish should be released in the downstream reach (0.5L < x < L) if the upstream part is not improved.