Patricio Farrell
email: patricio.farrell@wias-berlin.de

Julien Moatti
email: julien.moatti@inria.fr

Michael O'donovan
email: michael.odonovan@wias-berlin.de

Stefan Schulz
email: stefan.schulz@tyndall.ie

Thomas Koprucki
email: thomas.koprucki@wias-berlin.de

Keywords: Scharfetter-Gummel scheme, (In, Ga)N devices, numerical simulations, thermodynamic consistency

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Importance of satisfying thermodynamic consistency in optoelectronic device simulations for high carrier densities

Importance of satisfying thermodynamic consistency in optoelectronic device simulations for high carrier densities

1 Introduction

In recent years drift-diffusion simulations have offered a numerically attainable method for studying carrier dynamics in a wide variety of devices including light-emitting diodes (LEDs) [START_REF] Römer | Signature of the ideality factor in IIInitride multi quantum well light emitting diodes[END_REF][START_REF] Li | Localization landscape theory of disorder in semiconductors. III. Application to carrier transport and recombination in light emitting diodes[END_REF][START_REF] O'donovan | Impact of random alloy fluctuations on the carrier distribution in multi-color[END_REF], transistors [START_REF] Szymaǹski | 2-D Drift-Diffusion Simulation of Organic Electrochemical Transistors[END_REF][START_REF] Darwish | A drift-diffusion simulation model for organic field effect transistors: on the importance of the Gaussian density of states and traps[END_REF] and solar cells [START_REF] Ren | Exploring the Way To Approach the Efficiency Limit of Perovskite Solar Cells by Drift-Diffusion Model[END_REF][START_REF] Tress | Optimum mobility, contact properties, and open-circuit voltage of organic solar cells: A drift-diffusion simulation study[END_REF].

The physical interpretation of the model is quite straightforward: carriers in a device tend to diffuse from regions of high carrier density to low, and have drift motion due to an applied force such as an electric field in the device.

On the other hand, the numerical implementation of the model can carry pitfalls which lead to an incorrect description of the device behaviour as we will highlight below in detail. The purpose of this publication is therefore to show the practical unphysical implications of disobeying thermodynamic consistency. In the stationary case, thermodynamic consistency for discretized drift-diffusion equations can be defined by the demand that the zero bias solution coincides with the thermodynamic equilibrium. In the transient case, it is closely related to the fact that, for boundary conditions compatible with the thermodynamic equilibrium, the solution converges to this equilibrium when time tends to infinity. It is already known that disobeying this property causes non-physical dissipation in the steady state, see [START_REF] Bessemoulin-Chatard | A finite volume scheme for convectiondiffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme[END_REF]. However, the inconsistent discrete approximation of the numerical fluxes has also more direct consequences for the quasi Fermi potentials and is thus important for accurately describing the physics of a device in the frame of drift-diffusion simulations. If one inconsistently approximates the fluxes the quasi Fermi level will show a completely wrong behavior in for instance quantum well regions, which are at the heart of modern LEDs. This has a knock-on effect for the description of the carrier densities and thus also recombination and current-voltage (IV) curves, which we will illustrate with an (In,Ga)N quantum well structure, a material system of strong interest for energy efficient solid state lighting [START_REF] Humphreys | Solid-State Lighting[END_REF] and for which considerable effort has been undertaken to develop advanced carrier transport models [START_REF] O'donovan | From atomistic tight-binding theory to macroscale drift-diffusion: Multiscale modeling and numerical simulation of uni-polar charge transport in (In,Ga)N devices with random fluctuations[END_REF].

For the Boltzmann distribution, the classical Scharfetter-Gummel scheme [START_REF] Scharfetter | Large signal analysis of a silicon Read diode[END_REF] presents such a thermodynamically consistent scheme. Strictly monotonically increasing non-Boltzmann distribution functions lead to diffusion enhancement. Various extensions of the Scharfetter-Gummel scheme have been suggested to account for this effect, see [START_REF] Jüngel | Numerical approximation of a drift-diffusion model for semiconductors with nonlinear diffusion[END_REF][START_REF] Purbo | Numerical model for degenerate and heterostructure semiconductor devices[END_REF][START_REF] Stodtmann | Numerical simulation of organic semiconductor devices with high carrier densities[END_REF]. Unfortunately, these schemes are not thermodynamically consistent. In [START_REF] Koprucki | Discretization scheme for drift-diffusion equations with strong diffusion enhancement[END_REF][START_REF] Koprucki | Generalization of the Scharfetter-Gummel scheme[END_REF] a thermodynamically consistent generalization for Blakemore statistics (which is itself a special case of [START_REF] Eymard | A finite volume scheme for nonlinear parabolic equations derived from one-dimensional local dirichlet problems[END_REF]) is presented in the spirit of [START_REF] Scharfetter | Large signal analysis of a silicon Read diode[END_REF] by solving local Dirichlet problems. But this generalization requires solving local nonlinear equations during assembly and the iterative solution of the coupled system. It is therefore computationally prohibitively expensive. A computationally more affordable approach is presented in [START_REF] Chainais-Hillairet | A monotone numerical flux for quasilinear convection diffusion equation[END_REF]. On the other hand, in [START_REF] Bessemoulin-Chatard | A finite volume scheme for convectiondiffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme[END_REF] the author presents another extension of the Scharfetter-Gummel scheme using a proper average of the nonlinear diffusion guaranteeing thermodynamic consistency for a specific choice of the distribution function. An alternative interpretation of this approach based on averaging the diffusion enhancement for a very general class of statistical distribution functions was given in [START_REF] Koprucki | On Thermodynamic Consistency of a Scharfetter-Gummel Scheme Based on a Modified Thermal Voltage for Drift-Diffusion Equations with Diffusion Enhancement[END_REF]. Finally, the SEDAN scheme [START_REF] Abdel | Assessing the quality of the excess chemical potential flux scheme for degenerate semiconductor device simulation[END_REF][START_REF] Cancès | A numerical-analysis-focused comparison of several finite volume schemes for a unipolar degenerate drift-diffusion model[END_REF][START_REF] Yu | SEDAN III -A one-dimensional device simulator[END_REF] includes the nonlinearity in the drift instead of the diffusion part of the flux and thus also yields a thermodynamically consistent scheme.

The remainder of this paper is organized as follows: In Section 2, we describe the bipolar drift-diffusion model for charge transport in semiconductors. Its finite volume discretization including the flux discretizations is described in Section 3. The formal definition of discrete thermodynamic consistency is also presented in this section. We compare thermodynamically consistent and inconsistent schemes in Section 4 by studying the distribution of densities and quasi Fermi levels within an (In,Ga)N quantum well (QW) system, which is embedded in a p-i-n junction. Finally, we conclude in Section 5.

Drift-diffusion equations and diffusion enhancement

We briefly introduce a model based on nonlinear partial differential equations which describes bipolar charge transport in a semiconductor. More details can be found in [START_REF] Farrell | Drift-diffusion models[END_REF]. The dependence of the carrier densities n and p on the chemical potentials for electrons and holes η n and η p are described by a statistical distribution function F as well as conduction and valence band densities of states N c and N v via the state equations n = N c F(η n ) and p = N v F(η p ). Typical choices for the distribution function are F(η) = exp(η), the so-called Boltzmann approximation, or

F(η) = F 1/2 (η) = 2 √ π ∞ 0 E 1/2
e E-η +1 dE, namely the Fermi-Dirac integral of order 1/2 describing degenerate semiconductors.

The chemical potentials are related to the quasi-Fermi potentials of electrons and holes φ n and φ p via

η n = (q(ψ -φ n ) -E c )/(k B T ) and η p = (q(φ p -ψ) + E v )/(k B T ).
Here q denotes the elementary charge, ψ the electrostatic potential, k B the Boltzmann constant, T the temperature and E c and E v the conduction and valence band-edge energies. This model assumes that charge carriers behave as if they are in a bulk material, described by a 3-D density of states. In a slowly varying potential this is a valid description, however in a quantum well system the abrupt interface requires a more advanced treatment. To account for the quantum mechanical nature of electrons and holes one option is to solve the Schrödinger equation for the confining potential energy formed by V c,v = E c,v -qψ. This is a numerically demanding approach, as the Schrödinger equation is an eigenvalue problem which would need to be solved Thermodynamic consistency in optoelectronic simulations self-consistently coupled to Poisson and drift-diffusion equations. While such a calculation is feasible in 1D, extending this to 2D or 3D structures is numerically prohibitive. Therefore, in recent years significant efforts have been undertaken to establish methods that account for quantum corrections but are computationally cheaper [START_REF] Li | Localization landscape theory of disorder in semiconductors. III. Application to carrier transport and recombination in light emitting diodes[END_REF][START_REF] Ferry | The Effective Potential in Device Modeling: The Good, the Bad and the Ugly[END_REF]. One of these approaches is based on the so-called localization landscape theory (LLT), which allows to extract a (non-local) effective potential. It has been shown that LLT provides a good approximation of the single particle ground states, not only in square wells but also triangular wells which are relevant for systems with a polarization field -such as (In,Ga)N QWs [START_REF] Chaudhuri | Electronic structure of semiconductor nanostructures: A modified localization landscape theory[END_REF]. The resulting (effective) confining potentials, E eff c,v , exhibit band edges that are softened and approximate the finite extent of carrier wavefunctions1 . In the framework of LLT a linear system of equations is solved to determine the effective potential without introducing extra free parameters [START_REF] Filoche | Localization landscape theory of disorder in semiconductors. I. Theory and modeling[END_REF]. Our implementation of this method in conjunction with a drift-diffusion based carrier transport solver is discussed in more detail in [START_REF] O'donovan | Impact of random alloy fluctuations on the carrier distribution in multi-color[END_REF] and [START_REF] O'donovan | From atomistic tight-binding theory to macroscale drift-diffusion: Multiscale modeling and numerical simulation of uni-polar charge transport in (In,Ga)N devices with random fluctuations[END_REF]. In [START_REF] O'donovan | Impact of random alloy fluctuations on the carrier distribution in multi-color[END_REF] we have compared results of our quantum corrected drift-diffusion model employing LLT with the results of a commercially available software package utilizing a 'standard' Schrödinger-Poisson solver. Our findings indicate that LLT can produce results in good agreement with the fully coupled Schrödinger-Poisson-drift-diffusion solver, highlighting that LLT captures quantum mechanical corrections sufficiently. In recent years this method has been used to study transport behaviour of numerous semi-conductor structures include LEDs [START_REF] Lheureux | A 3D simulation comparison of carrier transport in green and blue cplane multi-quantum well nitride light emitting diodes[END_REF][START_REF] Lynsky | Barriers to carrier transport in multiple quantum well nitride-based c-plane green light emitting diodes[END_REF][START_REF] Lynsky | Improved Vertical Carrier Transport for Green III-Nitride LEDs Using (In,Ga)N Alloy Quantum Barriers[END_REF][START_REF] Römer | Luminescence distribution in the multiquantum well region of III-nitride light emitting diodes[END_REF], blocking layers [START_REF] Qwah K S, Monavarian | Theoretical and experimental investigations of vertical hole transport through unipolar AlGaN structures: Impacts of random alloy disorder[END_REF] and superlattices [START_REF] Tsai | Application of localization landscape theory and the k•p model for direct modeling of carrier transport in a type II superlattice InAs/InAsSb photoconductor system[END_REF]; moreover, these methods have been used successfully alongside experimental studies to gain insight into device behaviour [START_REF] Lynsky | Improved Vertical Carrier Transport for Green III-Nitride LEDs Using (In,Ga)N Alloy Quantum Barriers[END_REF][START_REF] Römer | Luminescence distribution in the multiquantum well region of III-nitride light emitting diodes[END_REF].

We model a bipolar semiconductor device as a domain Ω ⊂ R d where the carrier transport in a self-consistent electrical field is described by a system of partial differential equations. In the steady-state case this drift-diffusion system consists of Poisson's equation for ψ and continuity equations for electrons and holes:

-∇ • (ε 0 ε r ∇ψ) = q(C + p -n), ⃗ x ∈ Ω, (1)
-∇ • j n = -qR, ∇ • j p = -qR, ⃗ x ∈ Ω. (2)
Here, ε r is the relative permittivity, C is the net doping profile, and R = R(n, p) describes carrier recombination. Electron and hole current densities can be expressed in terms of quasi-Fermi potentials by

j n = -qµ n n∇φ n , j p = -qµ p p∇φ p , (3)
or for any strictly monotonic Fermi-like distribution function F(η) in driftdiffusion form

j n = qµ n U T g n N c ∇n -n∇ ψ - E c q and j p = -qµ p U T g p N v ∇p + p∇ ψ + E v q , (4)
where µ n and µ p denote the electron and hole mobilities, respectively, and U T = k B T /q is the thermal voltage. The factor g can be defined in terms of densities, g(x) = x(F -1 ) ′ (x), for x ∈ R. This factor is the so-called diffusion enhancement appearing as a density-dependent modification factor in the generalized Einstein relation, see [START_REF] Van Mensfoort | Effect of Gaussian disorder on the voltage dependence of the current density in sandwich-type devices based on organic semiconductors[END_REF], leading in general to a non-linear diffusion coefficient. For the Boltzmann distribution, F(η) = exp(η), we have g ≡ 1 and the current expressions (4) reduce to the usual ones with linear diffusion.

Finite volume space discretization

We discretize the domain Ω using the Voronoï box based finite volume method introduced in [START_REF] Macneal | An asymmetrical finite difference network[END_REF], also known as "box method" due to [START_REF] Bank | Some error estimates for the box method[END_REF]. It uses a simplical boundary conforming Delaunay grid ( [START_REF] Si | Boundary conforming Delaunay mesh generation[END_REF]) which allows to obtain control volumes surrounding each given collocation point x K by joining the circumcenters of the simplices containing it, see [START_REF] Farrell | Drift-diffusion models[END_REF] for details.

Let ∂K denote the boundary of the control volume K, and |ξ| the measure of a geometrical object ξ. For each control volume K, we integrate the continuity equation ( 2) and apply Gauss's theorem to the integral of the flux divergence. Restricting our considerations to the electron transport equation, we obtain

0 = ∂K j n • nds - K qR dx = L neighbor of K ∂K∩∂L j n • n KL ds - K qR dx ≈ L neighbor of K |∂K ∩ ∂L|j n,KL -q|K|R(n K , p K ), ( 5
)
where n is the internal unit normal to ∂K and n KL is the internal unit normal to the interface ∂K ∩ ∂L for each neighbor L of K. The values n K , p K are the numerical approximations of the densities n and p at the collocation points x K , and j n,KL are approximations of the normal currents through ∂K ∩∂L. In the same manner the discretization of the Poisson equation can be obtained. A more detailed discussion of this method can be found in [START_REF] Farrell | Drift-diffusion models[END_REF].

Discrete thermodynamic consistency

One property which holds on a continuous level to avoid unphysical state dissipation is the preservation of thermodynamic equilibrium [START_REF] Farrell | Drift-diffusion models[END_REF]. Mathematically, this means that vanishing fluxes shall imply constant quasi Fermi potentials. The classical discrete counterpart of this property is formulated as below (see for example [START_REF] Koprucki | On Thermodynamic Consistency of a Scharfetter-Gummel Scheme Based on a Modified Thermal Voltage for Drift-Diffusion Equations with Diffusion Enhancement[END_REF][START_REF] Farrell | Drift-diffusion models[END_REF]): a numerical flux j = j KL is said to be thermodynamically Thermodynamic consistency in optoelectronic simulations consistent if it satisfies an analogous discrete relation, i.e. j = 0 implies δφ KL = 0, [START_REF] Ren | Exploring the Way To Approach the Efficiency Limit of Perovskite Solar Cells by Drift-Diffusion Model[END_REF] where δφ KL = (φ L -φ K )/U T . Similarly, we define δη KL = η L -η K and δψ KL = (ψ L -ψ K )/U T and δE KL = (E c,L -E c,K )/(qU T ). We point out that the condition (6) holds in equilibrium. Here, we introduce a stronger notion of thermodynamic consistency, which holds outside of equilibrium, namely j ≤ 0 implies δφ KL ≥ 0 and j ≥ 0 implies δφ KL ≤ 0.

An important property of defining thermodynamical consistency like above is that the sign of the numerical current is consistent with that of its continuous counterpart [START_REF] O'donovan | Impact of random alloy fluctuations on the carrier distribution in multi-color[END_REF]. Thermodynamic consistency is also important, when coupling the van Roosbroeck system to heat transport models [START_REF] Farrell | Drift-diffusion models[END_REF]. We discuss now different numerical fluxes that may be used within a Voronoï finite volume framework.

The Scharfetter-Gummel scheme

First, we introduce the well known, classical Scharfetter-Gummel flux approximation [START_REF] Scharfetter | Large signal analysis of a silicon Read diode[END_REF] given by

j sg = j 0 {B (δψ KL -δE KL ) F(η L ) -B (-δψ KL + δE KL ) F(η K )} , (8)
where the constant j 0 is given by j

0 = qµ n N c U T h KL for h KL = |⃗ x K -⃗ x L |, and
B is the Bernoulli function, B(x) = x exp(x)-1 .
It is important to point out that Scharfetter and Gummel introduced this numerical flux only in the Boltzmann regime, i.e. F = exp. In this case, the flux is thermodynamically consistent in the sense of [START_REF] Ren | Exploring the Way To Approach the Efficiency Limit of Perovskite Solar Cells by Drift-Diffusion Model[END_REF]. However, once we leave the Boltzmann regime, i.e. F ̸ = exp, and continue using [START_REF] Bessemoulin-Chatard | A finite volume scheme for convectiondiffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme[END_REF] this numerical flux will no longer be thermodynamically consistent.

SEDAN scheme

Next, we present the SEDAN scheme, which yields a thermodynamically consistent approach also for state equations which do not necessarily rely on the Boltzmann approximation. The earliest reference we could find for such a excess chemical potential scheme is the source code of the SEDAN III simulator [START_REF] Yu | SEDAN III -A one-dimensional device simulator[END_REF], which explains the reason we use this name. A numerical analysis focused comparison of this flux approximation is given in [START_REF] Cancès | A numerical-analysis-focused comparison of several finite volume schemes for a unipolar degenerate drift-diffusion model[END_REF] and simulation results are presented in [START_REF] Abdel | Assessing the quality of the excess chemical potential flux scheme for degenerate semiconductor device simulation[END_REF]. The scheme is motivated by rearranging the drift part to include the excess chemical potential, ν ex = ln F(η) -η, yielding Fig. 1 Conduction band edge (black) and quasi Fermi energy (red) at a bias of 3.3V (a) when using Boltzmann statistics, (b) when incorrectly using the Scharfetter-Gummel (SG) scheme with Fermi-Dirac (FD) statistics, and (c) when correctly using the SEDAN scheme with Fermi-Dirac statistics.

j sedan = j 0 {B (Q KL ) F(η L ) -B (-Q KL ) F(η K )} (9)
with

Q KL = δψ KL -δE KL + ν ex L -ν ex K = δψ KL -δE KL -δη KL + ln F(η L ) F(η K ) . (10)
Note that, using the definition of Q KL and the fact that e x B(x) = B(-x), one can reformulate the SEDAN flux as

j sedan = j 0 B (Q KL ) F(η L ) 1 -e δφ KL .
Therefore, it is easy to see that the SEDAN flux satisfies [START_REF] Tress | Optimum mobility, contact properties, and open-circuit voltage of organic solar cells: A drift-diffusion simulation study[END_REF], since both B(Q KL ) and F(η L ) are positive: it is a thermodynamicaly consistent numerical flux.

Note that when applying Boltzmann statistics ν ex = 0 and the SEDAN flux becomes equivalent to the Scharfetter-Gummel expression. Therefore, in the next section, when displaying results using Boltzmann statistics, we only show results from one numerical scheme.

Simulations

To illustrate the importance of using thermodynamically consistent flux approximations, we study a simple (In,Ga)N multi QW (MQW) system. In particular, we consider three QWs and the same set of parameters as in [START_REF] O'donovan | Impact of random alloy fluctuations on the carrier distribution in multi-color[END_REF]. For large negative values of η (η ≤ -2, which correspond roughly to densities below 14% of the effective density of states N c , thus a low carrier density regime in the conduction band of the wells) the Boltzmann approximation provides a good estimate of the Fermi-Dirac statistics. Therefore in certain cases Scharfetter-Gummel scheme can offer a good description of the drift-diffusion model (e.g. Ref. [START_REF] O'donovan | Multiscale simulations of uni-polar hole transport in (In,Ga)N quantum well systems[END_REF]).

On the other hand there are situations where Boltzmann statistics will not suffice and the importance of using a thermodynamically consistent scheme becomes apparent. Figure 1 (a) displays the conduction band edge and quasi Thermodynamic consistency in optoelectronic simulations Fermi energies of the MQW system when treated using Boltzmann statistics at a bias of 3.3V. Applying the Scharfetter-Gummel scheme [START_REF] Bessemoulin-Chatard | A finite volume scheme for convectiondiffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme[END_REF] to Fermi-Dirac statistics for a bias of 3.3V leads to humps in the quasi Fermi energy within each QW, see Figure 1 (b).

Recalling the condition [START_REF] Tress | Optimum mobility, contact properties, and open-circuit voltage of organic solar cells: A drift-diffusion simulation study[END_REF], for thermodynamically consistent schemes, the discrete gradient of the discrete quasi Fermi level should indicates the direction of the discrete electron flow (as it is in the continuous case, according to (3)). However in Figure 1 (b), one can see that, when the Scharfetter-Gummel scheme is applied to Fermi-Dirac statistics, resulting in a thermodynamically inconsistent scheme, the direction of electron flow is to the right outside the QW regions (e.g. between -5 nm and 3 nm) but to the left inside the QW regions (e.g. between 10 nm and 15 nm), which is not in accordance with [START_REF] Tress | Optimum mobility, contact properties, and open-circuit voltage of organic solar cells: A drift-diffusion simulation study[END_REF]. Moreover, as there is no generation of carriers in the system this change in direction of the electron flux is highly unphysical. These humps in the quasi Fermi level are not present if one uses the Scharfetter-Gummel scheme with the Boltzmann approximation, which is a thermodynamically consistent scheme (Figure 1 (a)). Similarly, using Fermi-Dirac statistics with the SEDAN scheme does not exhibit the unphysical humps in the Fermi level (Figure 1 (c)).

Another perspective to interpret thermodynamic inconsistency is to note the incorrect interplay between quasi Fermi energies and local current fluxes. We see in Figure 2 that the local numerical electron fluxes are positive and decrease monotonically across the QW regions. This is true for all three settings, that is, Boltzmann statistics in combination with the classical Scharfetter-Gummel scheme, as well as Fermi-Dirac statistics using both consistent and inconsistent numerical schemes. By our definition of a thermodynamically consistent scheme (7), a positive local numerical flux should imply a negative quasi Fermi potential discrete gradient. In fact, this is true ) Position Fig. 3 Electron density at a bias of 3.3V (a) when using Boltzmann statistics, (b) when incorrectly using the Scharfetter-Gummel (SG) scheme with Fermi-Dirac statistics, and (c) when correctly using the SEDAN scheme with Fermi-Dirac statistics.

for both consistent settings as one can see in Figure 1 (a) and (c). However, for the inconsistent case the derivative of the quasi Fermi energies can become positive inside the wells, Fig. 1 Previous studies of the numerical flux approximations have shown that a thermodynamically inconsistent scheme can result in the incorrect sign of the particle flux [START_REF] Koprucki | On Thermodynamic Consistency of a Scharfetter-Gummel Scheme Based on a Modified Thermal Voltage for Drift-Diffusion Equations with Diffusion Enhancement[END_REF]. This is also reflected by the fact that (6) holds only for consistent schemes, which guarantee the physically correct sign of the current also far from equilibrium.

The physical reason of why an inconsistent scheme produces humps within the QWs becomes apparent when looking at the corresponding densities at a bias of 3.3V. The Fermi-Dirac function grows like a polynomial while the Boltzmann approximation grows for large densities exponentially, see for example Figure 50.9 in [START_REF] Farrell | Drift-diffusion models[END_REF]. This different behaviour leads to nonlinear diffusion, the so-called diffusion enhancement, for non-Boltzmann statistics, see (4), or [START_REF] Koprucki | On Thermodynamic Consistency of a Scharfetter-Gummel Scheme Based on a Modified Thermal Voltage for Drift-Diffusion Equations with Diffusion Enhancement[END_REF]. The Scharfetter-Gummel scheme neglects the diffusion enhancement assuming only linear diffusion, which has a knock-on impact on the carrier density: the densities calculated using Boltzmann statistics (Fig. 3 (a)) and using Fermi-Dirac statistics with the Scharfetter-Gummel scheme (Fig. 3 (b)) are visibly indistinguishable. However, in order to produce the same density between a Boltzmann and Fermi-Dirac calculation, the quasi Fermi levels must differ. This results in the unusual behaviour of the quasi Fermi level exhibited in Figure 1 (b). Comparing these densities with the correctly calculated density using Fermi-Dirac statistics in combination with the thermodynamically consistent SEDAN scheme (Figure 3 (c)) we see that the choice of statistics function will impact carrier density in the well, and more strongly in the barrier regions at the here chosen example bias of 3.3V.

Because it influences the carrier density, thermodynamic inconsistency has direct implications for the computed recombination rates as well as the currentvoltage (IV) curves. Next, we compare the recombination rates calculated using Boltzmann statistics with those calculated while incorrectly using the Scharfetter-Gummel scheme implementing Fermi-Dirac statistics. where R s i is the recombination rate associated with the process i (i ∈ {Shockley-Read-Hall, Radiative, Auger}) calculated with the scheme s (s ∈ {Boltzmann, SEDAN}). From this figure it becomes clear that the Boltzmann behaviour overall underestimates the recombination across the multi QW region of the device. In particular, the Auger recombination is underestimated by up to two orders of magnitude at a bias of 3.3V. These differences increase as the bias is increased (not shown). This can have consequences for overall device behaviour such as the internal quantum efficiency and the IV curves. The latter are shown in Figure 4 (b), where the decreased recombination current displayed in the Boltzmann and thermodynamically inconsistent Fermi-Dirac scheme leads to an underestimate of the current density by close to an order of magnitude at 3.6V.

The results highlighted above indicate that Fermi statistics implemented using a thermodynamically inconsistent scheme will result in Boltzmann-like behaviour in LED simulations -at least in terms of carrier and current densities. If this is extended to laser simulations the consequences can be even more dramatic, as the gain calculation depends on the difference between the electron and hole quasi Fermi energies [START_REF] Bandelow | Fabry-Perot Lasers: Thermodynamics-Based Modeling[END_REF], expressed by the so-called Fermi voltage. In this case the unphysical humps seen in Figure 1 (b) will lead to an incorrect prediction of the transparency density.

Conclusion

In this paper, we have shown the importance of using a thermodynamically consistent flux discretization when describing drift-diffusion processes within quantum well devices.

Using the classical Scharfetter-Gummel scheme with Fermi-Dirac statistics is an example of such an inconsistent scheme. Here we studied an (In,Ga)N multi quantum well structure as an example since it is a very important material system for optoelectronic devices. In this case, the Fermi levels show humps within the quantum wells resulting in an unphysical description of the direction of the current, e.g. assuming the usual continuous expression. This is explained by the omission of diffusion enhancement from the numerical current expression, that leads to a similar density distribution as using Boltzmann statistics. This has a knock-on effect for recombination and current-voltage behaviour, where using Fermi-Dirac statistics with a thermodynamically inconsistent scheme may incorrectly predict a Boltzmann-like behaviour.

Contrarily, for a thermodynamically consistent scheme, such as the SEDAN scheme, these unphysical humps in the Fermi levels disappear and accurate current curves and recombination processes are predicted. Thus, thermodynamically consistent schemes are essential to address open questions, such as the efficiency drop in modern light emitting devices and to reliably guide their design.

Fig. 2

 2 Fig. 2 Numerical electron flux averaged over each atomic plane at a bias of 3.3V shown for Boltzmann statisitcs (black), Fermi-Dirac statistics incorrectly using the Scharfetter-Gummel (SG) flux discretization (red) and Fermi-Dirac statistics correctly using the SEDAN flux discretization (blue).

  (b), also in disagreement of how the flux should behave at the continuous level (see the relation between the flux and the quasi Fermi potential (3)).

Fig. 4

 4 Fig.4(a) Difference in magnitude of the Shockley-Read-Hall (SRH, black), radiative (red) and Auger (blue) recombination rates between Fermi-Dirac and Boltzmann statistics at a bias of 3.3V, calculated as described in the main text. (b) Current density-voltage (IV) curves using Boltzmann statistics (black, dashed), Fermi-Dirac statistics using the Scharfetter-Gummel scheme (SG, red) and Fermi-Dirac statistics using the SEDAN scheme (black, solid), shown on a log scale.

From here on within this work Ec,v shall refer to the band edge values which are modified to account for the effective potential, Ec,v ≡ E eff c,v .

Funding: This paper was supported by the Leibniz competition 2022 (UVSimTec, K415/2021), Science Foundation Ireland (Nos.17/CDA/4789 and 12/RC/2276 P2) and Labex CEMPI (ANR-11-LABX-0007-01).


G Chaput

É Prévost

Review of the Atlantic Salmon NAC Pre-Fishery Abundance (PFA) and Catch Advice Model for West Greenland

INTRODUCTION

Atlantic salmon (Salmo salar L.) from anadromous populations of eastern North America can undertake extensive feeding migrations generally throughout the North Atlantic. Salmon from many different stock areas can aggregate at certain periods near Greenland and in the Norwegian Sea and these feeding aggregations attracted the attention of fishermen and resulted in fisheries on mixed-stocks during the first and second years at sea. Sampling of the fishery catches, analyses of biological characteristics, and recaptures of salmon originally tagged in homewaters indicated that the salmon in these fisheries originated from a large number of salmon producing areas in the North Atlantic, from both eastern North America and the northeast Atlantic areas [START_REF] Saunders | Canadian Atlantic salmon recaptured near Greenland[END_REF][START_REF] Paloheimo | Reduction of Atlantic Salmon (Salmo salar) Catches in Canada Attributed to the Greenland Fishery[END_REF][START_REF] Reddin | Distribution and biological characteristics of Atlantic salmon (Salmo salar) at Greenland based on the analysis of historical tag recoveries[END_REF]. Recent genetic stock identification analyses of salmon of North American origin captured at West Greenland indicate the diverse regional contribution of salmon to these fisheries, including major (> 20%) contributions from three regional groups (Labrador, Quebec, Gulf) in eastern North America and minor (< 5% per) contributions from the three other regions (Bradbury et al. 2016b; Figure 1). The harvests were also predominantly of salmon in their second year at sea that had not previously spawned, and therefore represented recruitment to the spawning stock with the expectation that many of those that survived would have returned to homewaters as two-seawinter salmon (ICES 2016).

There was sufficient concern regarding the impact of the Greenland mixed stock fishery on salmon populations in homewaters of both eastern North America [START_REF] Paloheimo | Reduction of Atlantic Salmon (Salmo salar) Catches in Canada Attributed to the Greenland Fishery[END_REF] and in the northeast Atlantic that a regional fisheries management organisation (RFMO), called the North Atlantic Salmon Conservation Organization (NASCO), was formed by international convention in 1983 for the purpose of negotiating fishery agreements for the high seas Atlantic salmon fisheries. Signatories to the NASCO convention included parties from all salmon producing countries in the North Atlantic as well as parties (Denmark on behalf of Greenland and the Faroe Islands) that fished salmon in the high seas.

Being strictly an RFMO, NASCO did not have the scientific structure to assess the status of salmon populations in the North Atlantic and consequently solicited scientific advice from an independent science body, the International Council for the Exploration of the Sea (ICES). ICES was formed in 1902, received full international status by international convention in 1964, and is an independent scientific body supported financially and professionally by countries in the North Atlantic, including North America. ICES, in its scientific mandate, was involved from the beginnings of NASCO in 1984, initially in the development of fishery monitoring programs and in the assessment of fishing practices at Greenland (ICES 1984).

The ICES Expert Group that addressed the questions from NASCO is called the Working Group on North Atlantic Salmon (WGNAS). This working group has convened annually since 1984, and twice annually in the earlier years. The reports prepared by ICES WGNAS contain a large amount of historical information on the characteristics of homewater and mixed stock fisheries and catches, on population specific indices of abundance including returns and spawning escapements, on indices of return rates of smolts and depending on the year, information on specific topics of interest related to emerging threats or innovations to Atlantic salmon conservation and management. The WGNAS reports have grown in volume over time as the questions from NASCO increased but also because the WGNAS reports are unique compilations of the historical and contemporary information on fisheries, status and trends of Atlantic salmon in the North Atlantic. Some of the information used by ICES to develop summaries of stock status [START_REF] Chaput | Overview of the status of Atlantic salmon (Salmo salar) in the North Atlantic and trends in marine mortality[END_REF], that describe the origin of salmon in mixed stock fisheries [START_REF] Reddin | A history of identification to continent of origin of Atlantic salmon (Salmo salar L.) at west Greenland[END_REF][START_REF] Reddin | Distribution and biological characteristics of Atlantic salmon (Salmo salar) at Greenland based on the analysis of historical tag recoveries[END_REF]Bradbury et al. 2016b), that describe approaches used to develop region and continent wide indices of abundance [START_REF] Potter | Estimating and forecasting pre-fishery abundance of Atlantic salmon (Salmo salar L.) in the Northeast Atlantic for the management of mixed-stock fisheries[END_REF]) and that document specific models for the provision of catch advice [START_REF] Chaput | Provision of catch advice taking account of non-stationarity in productivity of Atlantic salmon (Salmo salar) in the Northwest Atlantic[END_REF]) have been published in the primary literature.

The statistical treatment of data and the development of models have advanced from the initial approaches documented in [START_REF] Ices | [END_REF]. The changes in statistical modelling approaches represent a paradigm shift in the consideration of data inputs and the previous treatment of derived values as equivalent to observations. In particular, the estimation of abundance by sea age group, modelling and catch advice components of the mixed stock fishery question have changed from disjoint and independent analysis components to a coherent inference and forecast approach that integrates all aspects of this dynamic into a sequential and related process from spawners, through survival, harvests in fisheries to returns to homewaters [START_REF] Olmos | Evidence for spatial coherence in time trends of marine life history traits of Atlantic salmon in the North Atlantic[END_REF]ICES 2021a). The most recent catch advice for the West Greenland fishery from ICES was provided for the fishing years 2021 to 2023 using the PFA forecast model that considered only the 2SW returns and 2SW spawners to six regions of NAC (ICES 2021b). The proposed version going forward, and which is reviewed in the ICES Benchmark process, is the life-cycle model (LCM) approach described by [START_REF] Massiot-Granier | Embedding stock assessment within anintegrated hierarchical Bayesian life cycle modelling framework: an application to Atlantic salmon in the Northeast Atlantic[END_REF], [START_REF] Olmos | Evidence for spatial coherence in time trends of marine life history traits of Atlantic salmon in the North Atlantic[END_REF]ICES (2021a); the LCM has yet to be implemented by ICES for inference, forecasting and catch advice.

PROGRESSION OF NAC PFA AND FORECAST MODELS USED AT ICES

For a number of years, ICES (e.g. 1994ICES (e.g. , 2015) ) provided catch advice for the West Greenland fishery using various models at the North American stock complex level to forecast PFA for the upcoming fishery year at Greenland based on the run-reconstruction approach and one or two continent-wide explanatory variables (Table 1).

A stock independent explanatory variable for PFA was proposed based on an association between catch rates of salmon at sea and sea surface temperature (SST). Using a demonstrated curvilinear relationship between catch rates (number of salmon captured per length of net and soak time) of salmon and SST, the habitat area favourable to salmon production was calculated based on SST values in an area north of 41ºN latitude and west of 21ºW longitude, an area encompassing 5.36 million km 2, that was considered to include most of the known locations at sea for salmon of North American origin. The relative index of the area suitable for salmon was developed by weighting the area within 2 deg squares at each temperature group by the catch rate for the same temperature group. Correlations between PFA values derived from run-reconstruction and time periods for the habitat index highlighted the winter period (January to March) in the PFA year as being the most strongly associated with reconstructed PFA.

For those initial years, absolute changes in PFA were related to absolute changes in the environmental variable, rather than the conventional approach of considering that the environment mediates survival (Table 1). This was later modified with models that considered both environmental and spawner indices as explanatory variables for PFA abundance and used the appropriate multiplicative error structure (ICES 2001) that takes into consideration that habitat acts on PFA through survival rather than on absolute abundance. A linear form of the model fitted the natural log of PFA relative to the natural log of spawners and habitat variables (Table 1). The retention of the habitat variable in the models that included a spawning stock explanatory variable was justified on the basis of studies showing synchrony over wide geographic areas of survival rates and that the winter period appeared to be a critical stage for post-smolt survival and maturation (Scarnecchia 1989;[START_REF] Reddin | Sea-surface temperature and distribution of Atlantic salmon in the Northwest Atlantic Ocean[END_REF][START_REF] Rago | A continental run reconstruction model for the non-maturing component of North American Atlantic salmon: analysis of fisheries in Greenland and Newfoundland-Labrador[END_REF]Friedland et al. 1998). The lagged spawner variable used was an incomplete index of the spawning stock, calculated as the sum of spawners from four of the six geographic areas of eastern North America and excluded the contribution of estimated spawners from a large area of eastern North America (Gulf) for the reason that inclusion of the entire spawner component of North America weakened the association between PFA and lagged spawners (ICES 2001).

ICES also reported on evidence of decline in productivity of salmon in the marine part of the life cycle as shown by declines in return rates of smolts to anadromous adults at a number of index sites in the North Atlantic and sustained declines in returns to homewaters despite important reductions in marine fisheries. As well, in the continent wide models that incorporated the habitat and lagged spawner indices as explanatory variables, the relative proportion of the explained variance in PFA due to the habitat variable was declining relative to that of the lagged spawner variable. In the 2001 assessment using the PFA time series from 1978 to 1999, the lagged spawner variable alone accounted for 79% of the variance of the log transformed PFA values whereas the habitat variable alone only accounted for 15% the variance of the log transformed PFA values (ICES 2001). In the following year's assessment with data from 1978 to 2000, the habitat variable alone accounted for only 11% of the variance of the log transformed PFA values with an increasing association between PFA and lagged spawners (ICES 2002).

In the 2003 assessment, with the closure of the commercial fisheries in Labrador in 1999, alternate estimates of spawners for Labrador had not been developed and the lagged spawner index that included Labrador was not available. An alternate lagged spawner index that included the five other regions was developed and used in the regression model for PFA. When this new index was used in the regression model, the overall model that included the habitat index was not significant (ICES 2003). However, an analysis of the sequence of PFA and the revised lagged spawner index revealed a temporal structure within the data set that had not appeared previously and that could not be accounted for by the model used in previous years. Two states of Atlantic salmon production become evident with a transition state from 1988[START_REF] Jensen | Atlantic salmon at Greenland[END_REF](ICES 2003;;[START_REF] Chaput | Provision of catch advice taking account of non-stationarity in productivity of Atlantic salmon (Salmo salar) in the Northwest Atlantic[END_REF]. To capture this dynamic, a model that incorporated a break into two time periods, termed phases, was fitted to the data. The position of the change between the high production phase and the lower more recent production phase was initially set to 1989 as this PFA year was perceived as the midpoint in the slide from a low spawner index and high PFA abundance to a high spawner index and unchanged PFA abundance (ICES 2003;[START_REF] Chaput | Provision of catch advice taking account of non-stationarity in productivity of Atlantic salmon (Salmo salar) in the Northwest Atlantic[END_REF]. When the PFA, lagged spawner and habitat index variables were fit with a phase shift model, the habitat variable was not significant but the lagged spawner variable with an intercept shift in the mean PFA levels was highly significant, accounting for 82% of the variance in log transformed PFA. However, to be used in a forecast context, an indication of the phase of productivity for the forecast year was needed. ICES (2003) used an independent approach outside the lagged spawer to PFA model to make a forecast of the productivity phase. Specifically, ICES indicated that it seemed reasonable to expect that productivity for the 2003 PFA would be in the lower phase, as observed over the last ten years. The approach taken to estimate this probability was to examine the historical changes in PFA from year t to year t+2, the two-year lag was used because current year PFA (i.e. 2002) was not available due to its dependence upon 2SW returns in the next year. Historical observations were used to estimate the possible values of PFA in the predicted year from the reconstructed PFA two years earlier to determine the probability of being in each phase of the two-phase regression. Consequently for the 2003 forecast of PFA, the probability of being in the high phase was estimated at 4.8% and the probability of being in the lower productivity phase was estimated as 95.2% (ICES 2003). The predicted PFA was a weighted combination of the two possible PFA distributions from two regressions, with weights determined by the probability of being in each phase. The two phase model for analyzing PFA and forecasting abundance was used by ICES for the subsequent years to 2010, including for multi-year advice in which the forecasts for 2009 and 2010 PFA years were based on the same models used for the 2008 PFA forecast with the probability of being in the low phase in 2009 and 2010 assumed to be the same as in 2008.

Alternate models that considered temporal changes in the mean productivity were also examined (ICES 2008) including a dynamic model, as described by [START_REF] Prévost | Static versus dynamic model for forecasting salmon pre-fishery abundance of the River Bush: a Bayesian comparison[END_REF]. There is no functional dynamic implied between PFA and LS in this model other than the temporal dynamic that production rate (𝛼 𝑦 synonymous with PFAy/LSy) in the year of interest would likely be similar to the previous year. The production rate (𝛼 𝑦 ) was modelled as either simple random walk or as an autocorrelated random shift between two productive states. In the exploration of these models, the mid-point values of PFA were modelled as a lognormal function of mid-point values of lagged spawners. The results confirmed the previous conclusions [START_REF] Chaput | Provision of catch advice taking account of non-stationarity in productivity of Atlantic salmon (Salmo salar) in the Northwest Atlantic[END_REF] in inferring that there had been a decrease in the productivity (PFA/LS) for the North American complex, beginning in 1989, continuing to decline into 1997, and remaining low thereafter (ICES 2008).

Forecasts of PFA for the fishery year of interest were derived using variants of models that always considered PFA at the continental scale as equivalent to an observation (Table 2). Some level of uncertainty regarding the relationship between PFA and the explanatory variables was eventually incorporated in the catch advice portion of the model by generating a number of data sets of estimated returns of 2SW salmon to North America (NR2.tot), catches of 2SW salmon in North America (NL2.tot) and catches of non-maturing 1SW salmon at Greenland (WG1nm.tot) and in North America (NL1nm.tot) based on random draws assuming uniform distributions of the input variables defined by minimum and maximum ranges from the run reconstruction. These simulated values of PFA, and eventually lagged spawners for the catch advice year of 2000, were used to summarize the distribution of the predicted PFA values from regression models adjusted to each stream of simulated data. Incorporation of measurement errors particularly in the predictor variable increased the uncertainty of the forecast. This is understandable because as the predictive variables become less informative about PFA, the most probable value of the PFA approaches the mean of the series. But this component of including uncertainty was not included implicitly in the inference component of the model.

Initially, the forecast of PFA was compared to the 2SW conservation requirements for eastern North America, adjusting for the losses due to natural mortality that would occur from the time of the PFA assessment to returns to homewaters (11 months of mortality). The risks of meeting the overall conservation requirement were thus provided using a risk analysis framework that considered the uncertainties in the PFA forecast and uncertainties in the characteristics of the fish in the fishery harvests. ICES (2000) recognized the inherent risks of managing a mixed stock fishery while not accounting for stock structure particularly when the status of the individual stocks differed greatly. Initial attempts to rectify this consisted of adjusting the overall spawning requirement to account for the risk of meeting conservation requirements simultaneously in the six regions of North America (ICES 2000;[START_REF] Chaput | Considerations for using spawner reference levels for managing single and mixed-stock fisheries of Atlantic salmon[END_REF]). This adjustment was the most optimistic scenario and assumed that rates of production were similar in all six stock areas. The reality was quite different as the stock status differed greatly among the regions with the expected returns of salmon to the USA and Scotia Fundy areas of recent years to be severely below their respective conservation requirements. The majority of the non-maturing 1SW salmon in the Northwest Atlantic in 2001 were expected to return principally to the Quebec, Gulf, Labrador and Newfoundland (ICES 2001). Ignoring these differences in anticipated relative production among regions and focusing on total abundance over regions did not provide a realistic assessment nor a useful analysis for guiding fisheries management.

In an attempt to rectify this, ICES (2002) partitioned the complex wide PFA forecast into potential regional returns based on the regional proportions of lagged spawners for the return year and assessed these regional returns against the region specific conservation requirements for 2SW salmon. Under estimated levels of 2SW spawning escapements, the majority of the non-maturing 1SW salmon in the Northwest Atlantic in 2001 were expected to return principally to the areas of Quebec, Gulf, Labrador and Newfoundland with very likely near zero chance for the Scotia-Fundy stock and zero chance for the USA stocks of meeting their region-specific 2SW spawner requirements. The partitioning of PFA to regions and returns to regions was done outside the stock dynamic model structure and in its approach assumed implicitly that there was one overall productivity dynamic common to all stocks from eastern North America.

Ultimately, ICES recognized the need to move to an alternate stock model that treated regions as specific components of the stock and recruitment dynamic which could be incorporated in a more cohesive and implicit structure with region-specificities in productivities and uncertainties. The model would implicitly include in the assessment of catch options the risk of achieving region specific spawner objectives.

NAC PFA MODEL

The onus to model region specific abundance relative to region specific spawning requirements initiated the reflection of incorporating implicitly the stock structure at the scale of regions within a single model. Also, it was recognized that PFA and lagged spawner inputs used in the complex wide models by ICES were not observations in the true sense, having been estimated from a number of other data sources, each with associated uncertainties, and generally these uncertainties were not included at the complex wide modelling scale (Table 2). ICES (2008) suggested that consideration be given to implicitly linking dynamics of multiple stages within a Bayesian modelling framework which can accommodate observation errors and parameter uncertainties. The regionally disaggregated model differs from the model used by ICES until then in the way observations are considered, the procedure for model fitting, and in the way inferences are drawn on the variables of interest. (Table 2) In the case of Atlantic salmon for the ICES context, the 2SW spawners and the 2SW returns are the closest inputs that could be considered observations. Intermediate observations include the catches which are modelled as covariates acting as controls on the abundance of salmon at different points at sea and ultimately on returns to homewaters. In the region-disaggregated model approach, the PFA (recruitment) is a latent variable. i.e. a variable conditioned by several components directly influencing its distribution (catches in the marine environment, returns to rivers, natural mortality) but which cannot be observed directly. The PFA stage is particularly problematic as there are no stage specific indices available which can be used to draw inferences on its value, and the observations that link to it are distant in time and space, either the spawners that produce it or the returns to homewaters that result from it. In the region-disaggregated model, it is the returns of 2SW salmon to regions that are treated as the observations, which in concert with fisheries catches at earlier stages at sea and assumptions on natural mortality post-PFA stage condition the estimates of abundance at the PFA stage. The region-disaggregated model structure is illustrated in Figure 2.

The recruitment rate dynamic between lagged spawners and returns is also modeled differently in the region disaggregated model from that previously used by ICES. The two phase model used by ICES (described by [START_REF] Chaput | Provision of catch advice taking account of non-stationarity in productivity of Atlantic salmon (Salmo salar) in the Northwest Atlantic[END_REF]) considers that there have been (and will be) two levels of recruitment rate experienced by the North American salmon complex. When the populations are in the low phase, they will either remain in the low phase or move to the high phase, there is no possibility of a further decline in recruitment rate or intermediate levels of recruitment rate. The region disaggregated model with a dynamic structure on the productivity parameter provides a more flexible structure that allows for increases or decreases over time regardless of the previous states of the populations (Prevost et al. 2005) although abrupt changes are not adequately detected because the annual changes are smoothed and the magnitude constrained by the relative changes estimated from the past.

NAC REGION DISAGGREGATED PFA MODEL

The NAC PFA forecast model used by ICES to 2021 builds on the previous work and expertise developed by the ICES North Atlantic Salmon Working Group (ICES 2009(ICES , 2012)). The purposes of the model are to :

• retrospectively estimate the NAC PFA prior to the marine fisheries on Atlantic salmon at Greenland and in the Newfoundland and Labrador sea fisheries based on "observed" lagged spawners, catches and returns to rivers (after marine fisheries), and

• to predict PFA for providing catch advice for the Greenland and Newfoundland and Labrador sea fisheries based on realized lagged spawners.

The PFA includes future 2SW fish only, referred to as 1SW non-maturing fish, and refers to the abundance prior to any sea fishery. It links the 2SW (lagged) spawner abundances in each of six regions to the subsequent 2SW returns to rivers in each of six regions, PFA being an intermediate and unknown state (latent variable). In the latest use of the model by the ICES WGNAS in March 2021, abundance prior to the fisheries is predicted for three years 2021 to 2023. These predicted abundances are used to provide catch advice for the fishery at West Greenland on 1SW non-maturing salmon in 2021 to 2023 and for the potential fisheries in North America on 2SW salmon when they return to homewaters for the years 2021 to 2024. The risk analysis and catch advice procedures have been described previously [START_REF] Chaput | Provision of catch advice taking account of non-stationarity in productivity of Atlantic salmon (Salmo salar) in the Northwest Atlantic[END_REF]ICES 2012) and are not presented here.

Materials and Methods

The North American stock complex is partitioned into six regions (k; Figure 1):

1. Labrador, 2. Newfoundland, 3. Quebec, 4. Gulf of St-Lawrence, 5. Scotia-Fundy,

USA

Outputs from the run reconstruction model for 2SW returns and 2SW lagged spawners [START_REF] Rago | Index measures and stock assessment in Atlantic salmon[END_REF]) by stock unit are used as inputs to the PFA forecast model (Figure 2).

The 2SW lagged spawers (LSi,k) represent the sum of smolt age adjusted annual 2SW spawners by stock unit (k) and year (i) that would be expected to contribute to the recruitment at sea prior to the fisheries (PFA) for year i (Figure 3). LSi,k are not directly observed but are estimated from the run-reconstruction submodel. The probability distributions of LS and returns of 2SW salmon by stock unit are used as likelihood functions expressing comparative degrees of belief given the data and a probability model not explicitly specified in the current model. The probability distributions were drawn from the run reconstruction Monte Carlo simulations and assumed normal with known mean (LS.m) and precision (1/variance) (tau.LS). The use of this distribution as a likelihood function is equivalent to assuming a pseudo-observation equal to LS.m issuing from a sampling distribution with mean and precision equal to LS and tau.LS [START_REF] Michielsens | Combining multiple Bayesian data analyses in a sequential framework for quantitative fisheries stock assessment[END_REF].

LS.mi,k ~ N (LSi,k, tau.LSi,k)

The LS.mi,k (mean) and tau.LSi,k (precision) were derived assuming the lagged spawner values were issued from a normal distribution characterized by the mean and standard deviation statistics from the run reconstruction Monte Carlo simulations of spawners and lagged spawners (Figure 4). For stock unit 6 (USA), the run reconstruction does not provide a dispersion for the Lagged spawners; the precisions (tau.LS[,6]) were drawn from a normal distribution assuming a fixed coefficient of variation of 0.01.

Returns of 2SW salmon back to the rivers of North America (NR2i,k) are not directly observed but estimated with uncertainty (Figure 5). The run reconstruction model provides a probability distribution of NR2 by region, taken as comparative degrees of credibility of NR2, conditional on observed data and expertise. The way this expertise is used to link NR2 to the observed data (e.g. river counts) is not included in the current model. The NR2 probability distributions from run reconstruction are used as likelihood functions assumed to be normal with known mean NR2.m and precision tau.NR2. The use of this distribution as a likelihood function is equivalent to assuming to have one pseudo-observation equal to NR2.m issuing from a sampling distribution with mean and precision equal to NR2 and tau.NR2 [START_REF] Michielsens | Combining multiple Bayesian data analyses in a sequential framework for quantitative fisheries stock assessment[END_REF]).

NR2.mi,k ~ N (NR2i,k, tau.NR2i,k)

The NRmi,k (mean) and tau.NR2i,k (1/variance) were derived assuming the returns were issued from a normal distribution characterize by the mean and standard deviation of the Monte Carlo simulations of returns. For stock unit 6 (USA), the run reconstruction does not provide a dispersion for the Lagged spawners or the 2SW returns. The precisions (tau.LS [,6] and tau.NR2[,6]) were drawn from a normal distribution assuming a fixed coefficient of variation of 0.01.

Modelling the lagged spawners to 2SW returns

The years are modeled independently conditionally on the lagged spawners and yearly productivity parameters (Figures 2,3). The lagged spawners to PFA ratios (productivity) are modeled dynamically, i.e. assuming they are sequentially dependent within a region and attempts to take into account the most significant sources of uncertainty. The DAG for the model is shown in Figure 3.

PFA is assumed to be proportional to lagged-spawners (LS), with i.i.d. lognormal errors, and is modeled separately for each region (k = 6). The first year in the time series (t) is 1978 for lagged spawners (due to the range of smolt ages 1 to 6 for NAC and the start of the spawner time series in 1970) and the last year of lagged spawner data is for the 2023 PFA year. In ICES (2021b), the PFA was modeled for 1978 to 2019 (the last PFA year for which returns of 2SW salmon have been estimated back to rivers in 2020).

²) . , ( , ,

PFA PFA LogN PFA k i k i  = k i k i k i k i a LS PFA , , , , ) log( 
+ + = ²) . , 0 ( , PFA N iid k i  
The total PFA is calculated as the sum of the regional PFA's (k = 6). The proportion of the total PFA in each region is calculated directly as: p.PFAi,k = PFAi,k / PFA.toti A non-informative prior is assumed for

The proportionality coefficient (log)

k i a , between LSi,k and PFAi,k for each region, also called the productivity parameter, is modeled dynamically as a random walk with a year and stock unit residual variation ( k i,  ) assumed multivariate normal (MVN). The variance covariance matrix (  ) allows for correlations among stock unit productivity values reflecting that the fish share a common marine environment during part of their life cycle and that there are regional specificities in the evolution of the freshwater or the marine coastal environment.

k i k i k i a a , , , 1  + = + ) , 0 ( ,  MVN iid k i



The common yearly evolution of ai is the mean of annual a across regions: a.yi <-mean (ai,k) This parameterization of the covariance of the proportionality coefficient is similar to the version of the model used by ICES (2015ICES ( , 2018)). The correlation matrix of a among the regions is calculated from the covariance matrix.

PFA to returns of 2SW fish to regions of North America

The catches in the commercial fisheries of West Greenland and the Newfoundland and Labrador commercial and FSC fisheries (NG1.tot, NC1.tot and NC2.tot) are not directly observed but estimated with error. The catches are converted to numbers of 1SW non-maturing and 2SW fish based on biological characteristics of the fish in the catch (Table 3). Their (prior) probability distributions are obtained from catch statistics according to a formal structure included in the current model.

The West Greenland fishery has two components: a reported catch and an unreported catch: For the Labrador and some of the Newfoundland fisheries (Table 4), 1SW salmon are present in the large salmon category, a portion of which are considered to be 1SW non-maturing. As well, it is assumed that some 1SW salmon in the small salmon catch category are also 1SW non-maturing.

WGHarv.Esti = WGHarvi + WGUnHarvi
For the fisheries in Labrador and SFA 3 to 7, the catch of 1SW non-maturing was estimated as: NC1.toti = f_imm * (LB_Com_Smi + NF3to7Com_Smi + (LB_Com_Lgi + NF3to7Com_Lgi) *q) + af_imm*(Nsm_LBFSCi + NLg_LBFSCi * q) with f_imm ~ uniform (0.1 to 0.2), af_imm ~ uniform (0.05 to 0.1), q ~ uniform (0.1, 0.3).

For the fisheries in Newfoundland and Labrador (SFA 1 to 14B), the catch of 2SW salmon was estimated as:

NC2.toti = (LB_Com_Lgi + NF3to7Com_Lgi + NLg_LBFSCi) * (1 -q) + NF8to14A_Lgi with q as above ~ uniform (0.1, 0.3) Escapement (N1i,k) from the Greenland (NG1.toti) and NFLD/Labrador (NC1.toti) fisheries on nonmaturing 1SW salmon are calculated by substracting catches of non-maturing 1SW fish from the regional PFAs apportioned according to p.PFAi,k (i.e. each region is exploited at the same rate).

N1i,k = PFAi,k -(NG1.toti + NC1.toti) * p.PFAi,k Exploitation rates of the regional stocks for any of the fisheries are estimated from the catch and the corresponding regional abundances at the time of the fisheries. Exploitation rate on 1SW non-maturing at West Greenland:

ER.WGi,k = NG1i,k / PFAi,k
Exploitation rate on 1SW non-maturing salmon in North America:

ER.C1i,k = NC1i,k / PFAi,k
Returns of 2SW fish (N2.1i,k) back to the coast of North America take into account natural mortality M1 over 10 months, the average delay between the Greenland and NFLD/Labrador fisheries on 1SW nonmaturing fish and the arrival of the 2SW fish along the coast of North America. M1 is assumed to be iid across years according to a normal distribution with mean 0.03 and 95% interval [0.02, 0.04].

N2.1i,k = N1i,k * si with si = exp -M1 * 10 and M1 ~ N(0.03, 2.5*10 -5 )

Without this narrowing of the potential range of variation of M1, there is a confusion between PFA and M1. This is a very strong a priori hypothesis that assumes no temporal variation in the average M1 such that all the temporal variation contained in the LS to 2SW returns must be accounted by the evolution over time of the  parameter. In addition, the smaller the value of  (a priori), the stronger the influence (a priori) of NC1 and NG1 on the returns to North America.

Regional apportioning of the 2SW catches in Newfoundland, Labrador and St. Pierre and Miquelon fisheries

Based on expert opinion, the exploitation rates of 2SW of the regional groups in North America differed within the fisheries of Newfoundland and Labrador. Specifically, the commercial and FSC fisheries of Labrador were assumed to exploit a higher proportion of Labrador origin 2SW salmon than would be the case based on the relative abundance of PFA to that time and those fisheries [START_REF] Bradbury | Genetic evidence of local exploitation of Atlantic salmon in a coastal subsistence fishery in the Northwest Atlantic[END_REF]. The following assumptions were made to apportion the catches of 2SW salmon in these fisheries to each region.

The fisheries of Labrador (SFAs 1, 2, 14B; Table 4) were assumed to be comprised of the following proportions:

• pLabi ~ uniform (0.6, 0.8) for i = 1971 to 1997

• pLabi ~ uniform (0.9, 1.0) for i = 1998 to present

• the non-Labrador portions of the catches were apportioned as pPFA'k for the other five regions (pPFA'i,k = pPFAi,k / sum(pPFAi,k, excluding Labrador)

The commercial fisheries harvests of Newfoundland (Table 4) were assumed to be comprised of the following regional proportions:

• the commercial fisheries of Newfoundland in SFA 3 to 7 (northeast coast) were apportioned to the regions as pPFAi,k (i.e. including Labrador)

• the commercial fisheries of Newfoundland in SFA 8 to 14A (south and west coast) were apportioned to the regions as pPFA'i,k (i.e. no Labrador origin fish in these fisheries)

The estimated catches of large salmon from the St. Pierre and Miquelon fisheries (SPMC2i; Table 3) were apportioned to the regions as pPFA'i,k (i.e. no Labrador origin fish in these fisheries) (Bradbury et al. 2016a).

Escapement from the North American coastal marine fisheries (N2.2i,k) is obtained by subtracting the 2SW catch (NC2.toti,k) from the 2SW returns (N2.1i,k).

N2.2i,k = N2.1i,k -NC2i,k -SPMC2i,k
Exploitation rates of the regional stocks for any of the fisheries are estimated from the catch and the corresponding regional abundances at the time of the fisheries.

Exploitation rate on 2SW salmon in North America (Newfoundland and Labrador) fisheries:

ER.C2i,k = NC2i,k / N2.1i,k
Returns of 2SW salmon (NR2) back to the rivers of North America take into account natural mortality (M2) over 1 month, the average delay between the NFLD/Labrador fishery and the entry in rivers for the 2SW fish.

NR2i,k = N2.2i,k * si with si = exp -M2 * 1 and M2 ~ N(0.03, 2.5*10 -5 )

FORECASTS AND CATCH OPTION SCENARIOS (ICES 2021)

ICES (2021b) provided a forecast for the year 2020 (maximum lagged spawner year (2023) -3) for which the lagged spawners, and the catches of 1SW fish in the Greenland and NFLD/Labrador fisheries have been observed but the catch in the NFLD/Labrador fishery and the river returns of 2SW have not yet been observed. This part of the model can be used to provide catch advice on 2SW salmon for the NFLD/Labrador fishery in the year 2021.

Modelling for the years Y-2 (2021) to Y (2023), for which only the lagged spawners have been observed, is used to forecast PFA abundance by region and to provide catch advice for the Greenland and the NFLD/Labrador fisheries. River returns are predicted conditionally on a catch option for the Greenland fishery in year i and for the NFLD/Labrador fishery in the year i and i+1. Greenland catch options are expressed in weight (t) and a sharing arrangement for North America and West Greenland of 60:40 (effectively means that for any catch option examine for West Greenland, the catch used in the risk analysis is raised by 2.5 (WGcatch/0.4).

Risk analysis is based on the management objectives being simultaneously achieved in the six regions (ICES 2021b). The biological characteristics of the salmon in the fishery at West Greenland are estimated based on the characteristics of the most recent five years (Table 3). Biological characteristics for the fishery are drawn independently from a uniform distribution between the minimum and maximum values of the recent five years.

RETROSPECTIVE ANALYSIS OF MODEL

The forecast abundances for PFA and the productivity parameter for the forecast years 2018 to 2020 (from ICES 2018) were contrasted with the reconstructed values based on the 2SW returns to the six regions up to 2020. The comparison that was conducted by ICES (2021b) is presented in this section.

RESULTS AND DISCUSSION

Lagged spawners overall for NAC have generally been less than half the 2SW conservation limit for NAC (Figure 4). The lowest lagged spawner values were estimated during the 2003 to 2013 PFA years, with a slight improvement in abundance for the 2015 to 2016 and higher values for the 2020 to 2023 PFA years. The improvements in 2SW spawners in Labrador during 2013 to 2017 are now accounted for in the lagged spawners and these are the major contributors to the increased number of lagged spawners for NAC in the 2020 to 2023 PFA years.

Estimates returns of 2SW salmon to each region and overall for NAC for the years 1971 to 2020 (PFA years 1970 to 2019) are shown in Figure 5. Returns of 2SW salmon to rivers (after marine fisheries) exceeded the region-specific 2SW CLs at various times in the past, with exception to Labrador for which returns have only exceeded 2SW CL since 2012. The returns of 2SW salmon to US have never exceeded the 2SW CL for the region (29,990s fish) but exceeded the management objective for 2SW salmon in a few years in the 1980s (Figure 5).

PFA TO LS PRODUCTIVITY

The productivity coefficient (log of PFA to LS) was highest in most regions prior to 1990 (PFA year) and decreased in all regions to reach the lowest values in the mid 1990s and early 2000s (Figures 6 and7). Productivity coefficient values near zero or negative (a value of zero means the number of fish produced at the PFA stage equals the number of lagged spawners that produced them whereas a negative value means the PFA estimate was less than the lagged spawners) were estimated for Labrador and Newfoundland in the early 2000s, for Gulf during 1998 to 2000, and for Scotia-Fundy during the 1990s and again in the mid-2010s PFA year. The most recent year values (2019 PFA year) are positive for all regions (Figures 6 and7).

The productivity coefficient for NAC overall was negative in 2001 and improved from that point onward but remains at values less than half those estimated between 1978 and 1988 PFA years (Figures 7 and8).

There are weak correlations between the regional productivity parameters (Table 5).

PFA ESTIMATES

Estimated abundance at the PFA stage (August 1 of the second summer at sea) by region show similar temporal trends of higher abundance in the late 1970s through the 1980s, with a rapid decline in the late 1980s to early 1990s (Figure 9), with exception of Newfoundland for which PFA estimated abundance declined in the late 1990s to early 2000s. There is an increase in the PFA to 2014 for Labrador resulting directly from the increased 2SW lagged spawners and increasing productivity but the subsequent PFA abundance declined associated with a decline in the productivity (Figure 6). PFA declined and remains near historical lows in all other regions (Figure 9).

The regional contributions to the overall NAC PFA were relatively stable over the period 1980 to 2008 with over 70% of total PFA contributed by Quebec and Gulf regions, followed by Labrador with over 20% of the overall PFA (Figure 10). The Scotia-Fundy region contributed as much as 20% of the PFA for the 1984 PFA year but through the 2000s, has represented less than 5% of the total PFA and the US has never represented more than a few percentage of the total (Figure 10). The proportion of NAC PFA has increased substantially for Labrador since 2012 (Figure 10).

Forecasts of productivity and PFA for PFA years 2021 to 2023

The productivity parameter used in the forecast is the value derived from the last year in the model, with increasing uncertainty for each year of the forecast.

The overall productivity estimate for NAC in the most recent year PFA (2019) increased to a high positive value (median = 0.65; 1.9 fish at the PFA stage per lagged spawner) equal to that estimated previously during the 2007 to 2009 PFA years (Figure 8). By region, the most recent year value for the productivity was improved relative to the previous decade for Quebec, Gulf, Scotia-Fundy and USA while it remained low for Labrador or equal to values from the previous decade for Newfoundland. In all regions, the productivity value is low but positive compared to the estimates of the 1980s (Figures 6 and7).

For 2021 to 2023 PFA years, the 5th percentiles of the posterior distributions of the regional PFAs are less than the management objective reserves for Scotia-Fundy and USA. In addition, the 25th percentiles are below the objectives for Gulf (Figure 9). For NAC overall, the predicted values (25 th percentile) for 2021 to 2023 are all substantially below the 2SW CL reserve (CL corrected for 10 months of M = 0.03 per month; Figure 8). The forecasts have very high uncertainty and the uncertainties increase as the forecasts move farther forward in time.

RETROSPECTIVE ANALYSIS OF MODEL

In 2018, the ICES Working Group provided forecasts of the regional productivity parameters and the regional specific PFAs based on the regional lagged spawners. The productivity parameter used in the forecast is the value derived from the last year in the model, with increasing uncertainty for each year of the forecast. In the 2018 assessment, the productivity parameter used to forecast the 2018 to 2020 PFA years was negative for three regions (Gulf, Scotia-Fundy, USA), positive and at low values for Quebec and Newfoundland, and high for Labrador (ICES 2018; Figure 11). The returns of 2SW salmon in 2018 to 2020 were slightly higher than expected in all regions except Labrador and the realized productivity for the 2017 to 2019 PFA years was higher than predicted in 2018 (Figure 12). As a result the estimated regional PFA values were lower in Labrador for the 2017 to 2019 PFA years and slightly higher in all the other regions (Figure 13). The larger overestimate for Labrador relative to the other regions resulted in a lower PFA value for NAC for those years than forecast in the 2018 assessment. Due to the large uncertainty associated with the forecast values, the estimated PFA values for 2017 to 2019 were within the 95% confidence intervals of the forecast values.

The previous advice provided by ICES (2018) indicated that there were no mixed-stock fishery catch options on the 1SW non-maturing salmon component for the 2018 to 2020 PFA years and the 2021 assessment confirmed that advice. -Lagged spawners (sum of Newfoundland, Quebec, Gulf, Scotia-Fundy, USA, excluding Labrador) -Phase shift with eight candidate break years (1986 to 1993) -Incorporated uncertainty in PFA and LS variables as derived random draws from inputs to derive PFA and LS. -For each random data set, most parsimonious model was selected using Akaike's Information Criterion and the selected model was used to generate a value for the probability density for the forecast year of interest. -ICES (2004,2005,2006,2007,2008,2009) ²) . , ( , ,

TABLES

PFA PFA LogN PFA k i k i  = ; k i k i k i k i a LS PFA , , , , ) log( 
+ + = k i k i k i a a , , , 1  + = + ; ) , 0 ( ,  MVN iid k i



-Region disaggregated lagged spawners and 2SW returns -Bayesian model -ICES (2009,2012,2015,2018,2021) LS

[i,k] PFA[i,k] PFA.tot[i] a[i,k] σ.PFA a[i-1,k] [k] N1[i,k] NG1[i,k] Gld NC1[i,k] NFLB NG1.tot[i] Gld NC1.tot[i] NFLB N2[i,k] NC2[i,k] NFLB NC2.tot[i] NFLB S1[i] NR2[i,k] S2[i] Region k Year i Pseudo Likelihood NAC PFA Model Vers. 2012 SPMC2 [i] Pseudo Likelihood SPMC2 [i,k]

  The catch in weight (t) is converted into numbers of fish based on a sampled mean weight (kg): WGN.Esti = WGHarv.Esti * 1000 / WGMeanWti Only a portion of the catch at West Greenland is of North American origin (pNAC.est) and of nonmaturing 1SW salmon (p1SWNAC) : pNAC.esti ~ uniform (WGPropNACMini, WGPropNACMaxi) (for West Greenland fishery years 1978 to 1999, 2001; based on discriminant functions of continent of origin scale characteristics) pNAC.esti ~ beta (WGSampleNACi, WGSampleNEACi) (for West Greenland fishery years 2000, 2002 to present; based on genetic sampling) NG1.toti = WGN.Esti * pNAC.esti * p1SWNACi

FIGURESFigure 1 .

 1 FIGURES

Figure 2 .

 2 Figure 2. Summary of flow of data from run reconstruction to the NAC PFA inference and forecast model for2SW returns.

Figure 3 .

 3 Figure 3. Directed Acyclical Graph of the North American inference and forecast model for 2SW salmon.

Figure 4 .

 4 Figure 4. Median (2.5 th to 97.5 th percentile range) of spawners (circles) and lagged spawners (squares) of 2SW salmon to NAC overall and for each of the six regions. For spawners, year corresponds to the year of spawning. For lagged spawners, year corresponds to the year of PFA. The dashed line is the corresponding 2SW Conservation Limit for NAC overall and for each region; the 2SW CL for USA (29 990 fish) is off the scale in the plot for USA. The dotted line in the Scotia-Fundy and USA panels are the region-specific 2SW management objectives.

Figure 5 .

 5 Figure 5. Regional estimates of returns (grey circle) and spawners (open square) of 2SW salmon 1971 to 2020 (ICES 2021). The medians and 90% confidence interval ranges are shown. The dashed blue line is the corresponding 2SW conservation limit for each region. For Scotia-Fundy and US the dotted red line corresponds to the 2SW management objectives.

Figure 6 .

 6 Figure 6. Region specific PFA to LS ratio for PFA years 1978 to 2019. The values for 2020 to 2023 are predicted based on Lagged Spawners.

Figure 7 .

 7 Figure 7. Region specific PFA to LS ratio for PFA years 1978 to 2017. The values for 2020 to 2023 are predicted.

Figure 8 .

 8 Figure 8. Total PFA for NAC (top panel) and PFA to LS ratio (log; bottom panel) prior to exploitation. The dashed blue line in the top panel is the corresponding 2SW conservation limit corrected for 10 months of M between the PFA period and returns to homewaters for NAC.

Figure 9 .

 9 Figure 9. Region specific PFA values for PFA years 1978 to 2023. The values for 2020 to 2023 are predicted based on Lagged Spawners. The dashed blue line is the corresponding 2SW conservation limit corrected for 10 months of M between the PFA period and returns to homewaters for each region. For Scotia-Fundy and US the dotted red line corresponds to the 2SW management objectives corrected for M.

Figure 10 .

 10 Figure 10. Proportion of PFA in each region relative to overall PFA for NAC.

Figure 11 .

 11 Figure 11. Comparison of the estimated (median value) productivity parameter by region and overall for NAC (mean of regional values, black line) from the assessment in 2018 (upper panel; ICES 2018) and the corresponding productivity values estimated with updated values in the assessment from ICES (2021; (lower panel) for the PFA years 1978 to 2020. The points in both panels in the shaded rectangle are forecast values for the productivity parameter for the corresponding assessment periods.

Figure 12 .Figure 13 .

 1213 Figure12. Comparison of the estimated (median value) productivity parameter by region and overall for NAC (mean of regional values, black line) from the assessment in 2018 (x-axis; ICES 2018) and the corresponding productivity values estimated with data to 2020 (y-axis; ICES 2021). The symbols are colour coded for year (red early, blue later) and the forecast values fromICES (2018) and the estimated values from ICES 92021) are shown as green square symbols.

Table 1 .

 1 Progression of North American stock complex level models for PFA used by ICES. 𝑃𝐹𝐴 𝑦 ) = 𝛼 + 𝛽𝑃ℎ 𝑖 + (𝛾 + 𝛿 𝑃ℎ 𝑖 ) 𝑙𝑛(𝐿𝑆 𝑦 ) + 𝜀 Nested models including two without phase shifts and five models with phase shifts and with eight possible break year points (1986 to 1993) for each model.

	Model

NAC 2008 to present Input from run reconstruction

				NAC
				PFA C,a,t-i
		-Returns (2SW) by region (su)
	likelihood 2SW returns stock unit link in productivity	-Lagged spawners (2SW) by
		region (su)
					e -M
			Mar_C C,a,t-i
				Ret C,a,t
				 su
	Direct inputs to NAC PFA modelling and forecast		Ret su|C,a,t
	-Marine fisheries catches (NFLD, LAB, SPM, WG) -Fisheries catch characteristics (p1SWnonmat, …..)	HW_C su|C,a,t
	-M at sea -Time at sea between stages -CL 2SW spawner	Lag_Sp C,t+f+1  su	Lag_Sp su,t+f+1	Sp su|C,a,t Sp C,a,t  su	jurisdictions (su)
			pSmAge su,f

Forecast capacity

Forecasts are based on lagged spawner values available for three years beyond the last observed 2SW return year and an estimate of the likelihood of being in the high phase or the low phase of productivity. Forecast values take one of two levels of recruitment rate.

Same forecast capacity as Working

Group model excluding the need to estimate the probability of being in a high or low phase. Forecasts are based on estimated lagged spawners and the recruitment rate from the last observed year with variance from the entire time series.

Risk analysis of catch options

Assume characteristics of the catches will be similar to the range of values observed during previous five years. Catch options scenarios are explored..