BOATSv2: new ecological and economic features improve simulations of high seas catch and effort

Climate change and industrial fishing are having profound effects on marine ecosystems. Numerical models of fish communities and their interaction with fishing can help assess the biogeochemical and socioeconomic dynamics of this coupled human–natural system and how it is changing. However, existing models have significant biases and do not include many processes known to be relevant. Here we describe an updated version of the BiOeconomic mArine Trophic Size-spectrum (BOATS) model for global fish and fishery studies. The model incorporates new ecological and economic features designed to ameliorate prior biases. Recent improvements include reduction of fish growth rates in iron-limited high-nutrient low-chlorophyll regions and the ability to simulate fishery management. Features added to BOATS here for the first time include (1) a separation of pelagic and demersal fish communities to provide an expanded representation of ecological diversity and (2) spatial variation of fishing costs and catchability for more realistic fishing effort dynamics. We also introduce a new set of observational diagnostics designed to evaluate the model beyond the boundary of large marine ecosystems (66 commonly adopted coastal ocean ecoregions). Following a multi-step parameter selection procedure, the updated BOATSv2 model shows comparable performance to the original model in coastal ecosystems, accurately simulating catch, biomass, and fishing effort, and markedly improves the representation of fisheries in the high seas, correcting for excessive high seas and deep-sea catches in the previous version. Improvements mainly stem from separating pelagic and demersal energy pathways, complemented by spatially variable catchability of pelagic fish and depth- and distance-dependent fishing costs. The updated model code is available for simulating both historical and future scenarios.

1 Introduction

Recent developments have enabled the formulation of size-based fish community models based on fundamental ecological principles (Heneghan et al., 2021). Instead of resolving linkages between species or functional groups within marine food webs, aggregated size spectrum models are based on properties that emerge at higher levels of organization. These models rely on macroecological principles to connect individual growth and metabolism (Brown et al., 2004; Kooijman, 2010; Hatton et al., 2021) to community-level production and biomass (Gascuel et al., 2011; Blanchard et al., 2012; Maury and Poggiale, 2013; Jennings and Collingridge, 2015; Petrik et al., 2019; Heneghan et al., 2020). By simplifying complex ecosystem dynamics into community-level biodiversity (Maury, 2010; Petrik et al., 2019) and regional variations in trophic efficiency (Du Pontavice et al., 2020) and other ecological variables, these models can project the response of global marine ecosystems to warming and shifts in primary production due to climate change (Lotze et al., 2019; Tittensor et al., 2021).

The BiOeconomic mArine Trophic Size-spectrum (BOATS) model is a size spectrum model that incorporates an explicit representation of commercial fishing effort (Carozza et al., 2016, 2017). The model's integration of ecological and economic dynamics enables a clear illustration of the profound effects of advances in fishing technology on historical changes in fish biomass compared to the impacts of climate change (Galbraith et al., 2017; Carozza et al., 2019). The ability to simulate how fish catches respond to dynamic fishing effort allows optimization of the model's ecological parameters against observational reconstructions of fish catches (Carozza et al., 2017). Based on this optimization method, BOATS provided new estimates of the global biomass and cycling rate of fish, indicating a non-negligible impact of fishing on carbon sequestration and biogeochemical cycles (Bianchi et al., 2021; Le Mézo et al., 2022). It also offered mechanistic insights into the historical progression of fisheries (Guiet et al., 2020). While the model was originally designed under the assumption of open-access fishing effort, subsequent developments enabled it to investigate the effects of regulatory measures on fish community dynamics and their response to long-term and abrupt climate perturbations (Scherrer and Galbraith, 2020; Scherrer et al., 2020).

These studies prove the usefulness of BOATS for exploring various aspects of global fisheries. Still, comparisons with observations have also revealed discrepancies that suggest limitations in the model's parameterizations and missing mechanisms. For instance, high-nutrient low-chlorophyll (HNLC) regions are characterized by relatively low primary production despite available macronutrients (Moore et al., 2013). These regions represent more than one-quarter of the open-ocean surface area and include the Southern Ocean, the eastern equatorial Pacific, and the subarctic North Pacific. In HNLC regions, comparison of simulated effort with global reconstructions suggested excessive fishing activity in BOATS, indirectly pointing to excessive biomass accumulation in the model (Galbraith et al., 2019). Similarly, while the model provides a realistic representation of coastal fisheries, catches in the high seas (here defined as the regions of the ocean beyond large marine ecosystems) appear to be much larger than recent observational reconstructions. Note that large marine ecosystems (LMEs) are 66 coastal ocean regions defined by ecological criteria (Sherman and Duda, 1999). Specifically, BOATS simulates 40 % of global catches beyond the boundary of LMEs by the 1990s. This is approximately 4 times larger than the value of 8 %–9 % from empirical estimates (Watson, 2017; Pauly et al., 2020). The large high seas catches coincide with excessive fishing in waters found above the deep seafloor. During the 1990s, the model's catch-weighted mean depth of waters where fishing occurs is 1698 m, contrasting significantly with the observational reconstruction range of 154–372 m. This discrepancy limits the model's applicability to study the interaction of industrial fishing with high seas and deep-ocean ecosystems and suggests potential shortcomings in the representation of open-ocean food webs.

In parallel, recent studies have shed new light on large-scale aspects of global marine ecosystems and fisheries. Reconstructions of industrial fishing effort by Global Fishing Watch (GFW; Kroodsma et al., 2018) highlighted spatial variations in fishing costs (Sala et al., 2018) and revealed the importance of seamounts in concentrating fishing activity, especially for pelagic fisheries in the high seas (Kerry et al., 2022). New reconstructions of fishing effort that include artisanal and industrial sectors provide more nuanced insights on the development of regional fisheries (Rousseau et al., 2019, 2024). Regional catch reconstructions have revealed the importance of “straddling” species, which forage both within exclusive economic zones and in the high seas over the course of a year, thus disconnecting regions of fish biomass production from regions of biomass extraction (Sumaila et al., 2015). Analysis of catch data shows how energy inputs at the base of food webs determine the dominance of pelagic vs. demersal communities across latitudes (van Denderen et al., 2018), suggesting different temperature sensitivities of growth for these groups (van Denderen et al., 2020). Finally, harmonization and in-depth analysis of fishery-independent trawl data have begun to reveal large-scale fish biomass patterns with unprecedented accuracy (Maureaud et al., 2024).

Here, motivated by discrepancies between observations and simulations with the original BOATS model (BOATSv1) and insights from recent large-scale studies, we revise the model formulation to improve its representation of high seas vs. coastal fisheries and of pelagic vs. demersal communities, leading to a significant model update: BOATSv2. The rest of the paper consists of four main sections. Section 2 summarizes the main principles and formulation of BOATSv1. Section 3 details previous model developments and new features introduced in BOATSv2. Section 4 describes a revised model optimization procedure. Section 5 justifies the selection of an ensemble of five optimal parameters, compares the old and new model versions (highlighting improvements in the representation of global fisheries in BOATSv2), and discusses insights from the new formulation.

2 BOATSv1

The philosophy of BOATS is to ensure global applicability while including sufficient ecological and economic complexity to represent realistic first-order fishery dynamics. The model is designed to include a relatively small number of parameters and to be computationally efficient, facilitating objective parameter optimization. It uses vertically averaged habitat characteristics on a two-dimensional spatial grid to simulate the variability of fish communities, from small regions to the global ocean. In the following, we provide a brief overview of key model principles. We refer the interested reader to Appendix A for all equations and to previous publications for detailed explanations (Carozza et al., 2016, 2017). A schematic that illustrates the model principles (adapted for BOATSv2 from Carozza et al., 2017) is shown in Fig. 1.

Figure 1 Schematic diagram of the main modules, components, and processes of BOATSv2. Environmental forcings, shown in the left panel (“pelagic” for BOATSv1; “pelagic” and “demersal” for BOATSv2), drive an ecological module that solves for the evolution in time of fish biomass as a function of fish size for multiple groups with different maximum size, shown in the central panel. These fish biomass spectra interact with the dynamic of fishing, controlled by an economic module and economic forcings, shown in the right panel. Economic forcings are spatially uniform in BOATSv1 but can be spatially variable in BOATSv2. Environmental forcings include the spatial distribution of resources at low trophic levels (Π_ψ or Π_β) and representative habitat temperatures (T₇₅ or T_bot). Fish biomass spectra for multiple groups emerge from the balance of environmentally controlled growth (Γ, linked with ξ_P or ξ_VB), recruitment (R), natural mortality (Λ), and fishing mortality (H). Economic forcings, which include spatially uniform ex-vessel prices (p) and spatially variable fishing costs (c) and catchability (q), influence the dynamic of fishing effort (E) for each fish group. Color shades of forcings illustrate spatial variations, from low (light) to high (dark) values. This figure is updated from the schematic for BOATSv1 in Carozza et al. (2017).

2.1 Ecological module

BOATSv1 simulates the evolution in time (t) of the biomass of commercial fish as a function of size, f_k (in $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}$, where g represents grams of wet biomass; see also Table 1 for a list of variables and parameters), and its propagation along a spectrum of size classes [m₀,m_u], where m is the biomass of an individual fish (in g), and m₀ and m_u are the minimum and maximum fish sizes represented by the model. To include a coarse representation of species diversity, BOATS simulates n_k different fish “groups” distinguished by their asymptotic mass $m_{\mathrm{\infty },k}<m_{\mathrm{u}}$, labeled by the subscript k (see an illustration of the small and large groups in the central panel of Fig. 1). The total biomass density, $B_k={\int }_{m_{\mathrm{0}}}^{m_{\mathrm{\infty },k}}{f}_k\mathrm{d}m$, is the sum of each group's biomass across individual size classes (in g m⁻²). The propagation of biomass in each size spectrum as a function of time t is described by the McKendrick von Foerster equation:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial }{\partial t}}{f}_k=-{\displaystyle \frac{\partial }{\partial m}}{\mathit{\gamma }}_{S,k}{f}_k+{\displaystyle \frac{{\mathit{\gamma }}_{S,k}{f}_k}{m}}-{\mathrm{\Lambda }}_k{f}_k-h_k.\end{array}$$

The first term in Eq. (1) represents the rate of change in time of the fish biomass spectrum for each group. The second term is the divergence of the growth flux, i.e., the transfer of biomass to increasing size as fish grow. The third term encapsulates the biomass accumulation due to the increase in individual size as fish grow. The fourth and fifth terms represent losses from natural mortality and catch, respectively.

This first-order partial differential equation in time and size requires both a boundary condition, here prescribed at the smallest size class m₀ and representing recruitment, and an initial condition at t=0, representing the initial biomass distribution for each group:

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}{\mathit{\gamma }}_{S,k}{f}_k=R_{\mathrm{P},k}\frac{R_{\mathrm{e},k}}{R_{\mathrm{P},k}+R_{\mathrm{e},k}}\mathrm{for}m=m_{\mathrm{0}}\text{(boundary condition)}\\ f_k=f_{k,m,\mathrm{0}}\mathrm{at}t=\mathrm{0}\text{(initial condition)}.\end{array}\right.\end{array}$$

In Eqs. (1)–(2) the size-dependent growth rate γ_S,k (in g s⁻¹; Eq. 3) controls the biomass propagation through size (white arrows in the central panel of Fig. 1), influenced by local water temperature (T in Kelvin) and primary production Π_ψ (in $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). The natural mortality rate Λ_k (in s⁻¹; Eq. 4) represents the biomass losses within each size class from predation and natural mortality (gray line in the central panel of Fig. 1). The sink term h_k (in $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) is the biomass harvest by fishing that couples the ecological module to the fishery dynamics module (see Sect. 2.2). Finally, R_e,k and R_P,k in Eq. (5) are respectively the biomass input potential at the recruitment size m₀ from egg production (R_e) by mature individuals (see yellow arrows in Fig. 1) and the maximum biomass input potential at the recruit size given the primary production (R_P). Both modulate the total recruitment input, γ_S,kf_k (in $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). Simulations start from an initial biomass distribution $f_{k,m,\mathrm{0}}$ that approximates an ocean in the absence of fishing (“pristine”), estimated from environmental conditions (Sect. 3.3).

In BOATS, the growth rate at a given size occurs either at the maximum physiological rate when food is not limiting (gray area in the central panel Fig. 1) or proportionally to primary production Π_ψ when food is limiting (green area in the central panel Fig. 1). Accordingly, the growth rate is proportional to the minimum of two quantities: (1) the energy provided by primary producers that reaches a given size class ξ_P,k, given trophic transfer across the food web, divided by the number of fish in that size class, and (2) the maximum production potential for a fish of that size based on an individual-level allometric growth rate that follows a von Bertalanffy formulation ξ_VB,k (in g s⁻¹):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mathit{\gamma }}_{S,k}& =(\mathrm{1}-{\mathrm{\Phi }}_k){\mathit{\xi }}_{I,k}=(\mathrm{1}-{\mathrm{\Phi }}_k)min\left({\mathit{\xi }}_{\mathrm{P},k},{\mathit{\xi }}_{\mathrm{VB},k}\right)\\ & =(\mathrm{1}-{\mathrm{\Phi }}_k)min\left({\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{C},k}\mathit{\pi }m}{f_k}},{\mathrm{Am}}^{\mathrm{b}}-k_{\mathrm{a}}m\right).\end{array}\end{array}$$

Here, the term (1−Φ_k) accounts for a reduction of the biomass allocated to somatic growth, with a fraction Φ_k allocated to the generation of reproductive material, i.e., egg production (Eq. 5). Thus, when food is limiting, individual fish will grow according to $\mathit{\pi }={\mathrm{\Pi }}_{\mathit{\psi }}/m_{\mathit{\psi }}(m/m_{\mathit{\psi }}{)}^{\mathit{\tau }-\mathrm{1}}$, which defines a spectrum of available energy from primary production as a function of size (in $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$). Here, τ is the trophic scaling, and m_ψ is a representative cell size for primary producers (i.e., phytoplankton) at the base of the food web. The trophic scaling parameter determines the efficiency of the propagation of production through the consumer size spectrum to increasingly larger sizes and higher trophic levels, following the framework of the metabolic theory of ecology Brown et al. (2004). The representative size m_ψ is determined from the empirical phytoplankton size structure model of Dunne et al. (2005) and depends on temperature (T in °C) and primary production Π_ψ. To ensure coexisting fish groups and because of the scarcity of data available to constrain resource allocation, primary production is equally partitioned across the groups, i.e., ${\mathit{\varphi }}_{\mathrm{C},k}=\mathrm{1}/n_k=\mathrm{1}/\mathrm{3}$. While this is a first-order assumption that allows realistic simulation of catches by group, it should be revised as new observational constraints become available. When food is in excess of what can be consumed by the standing fish biomass, fish grow as fast as physiologically possible, given an allometric scaling b, a temperature-dependent anabolism Am^b=A₀ a_A(T)m^b, and catabolism $k_{\mathrm{a}}m=A{\mathit{ϵ}}_{\mathrm{a}}{m}_{\mathrm{\infty },k}^{b-\mathrm{1}}m$, where A₀ is a growth constant (in g s⁻¹) and ϵ_a an activity fraction. This formulation is inspired by an empirical allometric framework following the model of Von Bertalanffy (1949), where growth is determined by food intake after assimilation and standard metabolism, discounted from energy used in activity and reproduction. This maximum growth is temperature-dependent based on the factor a_A(T) (T in Kelvin), which follows a van't Hoff–Arrhenius curve controlled by a growth activation energy ω_a,A (in eV). The fish mortality is independent of variations in the growth rate.

The natural mortality rate (in units of s⁻¹) depends on both individual and asymptotic mass and represents biomass losses due to predation to organisms both within and outside of the resolved community size spectrum, as well as other natural causes. The natural mortality rate is based on an empirical parameterization (Gislason et al., 2010; Charnov et al., 2013):

$$\begin{array}{}\text{(4)}& {\mathrm{\Lambda }}_k=e^{{\mathit{\zeta }}_{\mathrm{1}}}{\displaystyle \frac{A_{\mathrm{0}}}{\mathrm{3}}}{a}_{\mathit{\lambda }}\left(T\right)m^{-h}{m}_{\mathrm{\infty },k}^{h+b-\mathrm{1}},\end{array}$$

where h is an allometric scaling, and ζ₁ (in g s⁻¹) a mortality rate parameter. As in Gislason et al. (2010), the natural mortality rate is linked to growth by means of the constants A₀ and b. To account for separate temperature dependencies between metabolism a_A(T) and other processes such as predator–prey interactions, mortality varies with a distinct temperature dependence a_λ(T), following a van't Hoff–Arrhenius curve controlled by a second activation energy ω_a,λ (in eV).

Recruitment provides the boundary condition, setting the flux of biomass at the lower mass boundary m₀. Recruitment is a function of the biomass production by mature individuals R_e,k and a limit to the survival of recruits that depends on resource availability, proportional to primary production R_P,k:

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{l}{R}_{\mathrm{e},k}={\mathit{\varphi }}_{\mathrm{f}}{s}_{\mathrm{e}}\frac{m_{\mathrm{0}}}{m_{\mathrm{e}}}{\int }_{m_{\mathrm{0}}}^{m_{\mathrm{\infty },k}}{\mathit{\gamma }}_{R,k}\left(m\right)\frac{f_k\left(m\right)}{m}\mathrm{d}m\\ R_{\mathrm{P},k}={\mathit{\varphi }}_{\mathrm{C},k}\mathit{\pi }\left(m_{\mathrm{0}}\right)m_{\mathrm{0}}.\end{array}\right.\end{array}$$

Assuming that a fraction Φ_k of the input energy is allocated to reproduction, ${\mathit{\gamma }}_{R,k}={\mathrm{\Phi }}_k{\mathit{\xi }}_{I,k}$, the recruitment flux is determined by summing the contributions of all mature individuals across sizes m for a fraction of females ϕ_f, an egg survival probability s_e, and a mean egg mass m_e (in g). Survival of recruits saturates towards a limit set by the energy available from primary production at the recruit size m₀. In high-biomass regions with large egg production rates, recruitment is thus generally constrained by π(m₀).

With this formulation, for a given set of parameters (see the list in Table 1 and more details in Appendix A), the dynamics of commercial fish biomass in BOATSv1 are determined by two spatially and temporally varying environmental “forcings” (shown as “pelagic” forcings Fig. 1): local primary production Π_ψ and epipelagic temperature T=T₇₅ (average temperature in the top 75 m in °C or K).

2.2 Economic module

BOATSv1 couples fish biomass and fishing effort E_k (in W m⁻²) to determine fish catch rates C_k (in $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) for each species group k (see harvest in the central panel of Fig. 1). Fishing effort is typically initialized everywhere at negligible values, starting from an unfished ocean, and evolves independently in each grid cell under an open-access dynamic, proportional to the local net profit (the difference between revenue and cost) as

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}t}{E}_k={\mathit{\kappa }}_{\mathrm{e}}\frac{{\mathrm{revenue}}_k-{\mathrm{cost}}_k}{E_k}\\ E_k=\mathrm{0}\mathrm{at}t=\mathrm{0},\end{array}\right.\end{array}$$

where κ_e (in ${\mathrm{W}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{USD}}^{-\mathrm{1}}$) is a fleet dynamic parameter that sets the effort adjustment timescale for a given change in profit. This formulation assumes an absence of regulation so that fishers seek the greatest total catch in each grid cell. When profitable, revenues in Eq. (7) exceed costs in Eq. (9), and fishing develops continuously. In the presence of a continuous increase in catchability, this eventually leads to a peak in catch, overharvesting, and ultimately biomass collapse (see Carozza et al., 2017, and Galbraith et al., 2017, for details).

The rate of revenue for a time interval dt from a given location in the ocean (in $\mathrm{USD}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is determined as a spectrum r_k (in $\mathrm{USD}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{g}}^{-\mathrm{1}}$) that is defined separately for each species group f_k. This represents the product of catch and the mass-specific price fishers are paid at port, integrated across size classes within each group:

$$\begin{array}{}\text{(7)}& {\mathrm{revenue}}_k=\underset{m_{\mathrm{0}}}{\overset{m_{\mathrm{\infty },k}}{\int }}{r}_k\mathrm{d}m\mathrm{d}t=p_k{q}_k{E}_k\mathrm{d}t\underset{m_{\mathrm{0}}}{\overset{m_{\mathrm{\infty },k}}{\int }}{\mathit{\sigma }}_k{f}_k\mathrm{d}m,\end{array}$$

where p_k is the ex-vessel price (in USD g⁻¹) for each group. The catchability q_k (in ${\mathrm{m}}^{\mathrm{2}}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) per unit of nominal fishing effort encapsulates the ability of fishing effort to extract fish biomass from the ocean. This quantity reflects the inherent characteristics of the fish group, as well as the fishing technology adopted (including gear, navigation instruments, sonars, and communication) and accrued knowledge (Palomares and Pauly, 2019). This formulation depends on the fraction of the fish biomass σ_kf_k targeted by fishing, where σ_k is a size-dependent selectivity of the fishing gear used to target group k. The selectivity plays a fundamental role by distributing fishing effort across size classes. A variety of functional forms exist, and all avoid the smallest sizes. These can be generalized as either dome-shaped (e.g., gill nets, longlines) or sigmoidal (e.g., trawls, seines, or dredges). Here, we parameterize the selectivity as a sigmoidal curve around a target threshold mass $m_{\mathrm{\Theta },k}=d_{m_{\mathrm{\Theta },k}}{e}_{m_{\mathrm{\Theta },k}}{m}_{\mathit{\alpha },k}$, essentially reducing the fishing effort targeting the smallest size classes:

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_k={\left[\mathrm{1}+{\left({\displaystyle \frac{m}{m_{\mathrm{\Theta },k}}}\right)}^{-c_{\mathit{\sigma }}/\mathrm{3}}\right]}^{-\mathrm{1}},\end{array}$$

with c_σ a fishing selectivity slope. The target threshold mass is proportional to the maturity mass for each group m_α,k, with the parameter $e_{m_{\mathrm{\Theta },k}}$ accounting for uncertainty around this mass and $d_{m_{\mathrm{\Theta },k}}$ set to select mainly mature individuals (i.e., $d_{m_{\mathrm{\Theta },k}}=\mathrm{1}$).

Net profits are determined by subtracting costs from revenues. Similar to revenue, the cost is expressed as the average expenditure rate per time over an area of the ocean (in $\mathrm{USD}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). In reality, the cost of fishing includes the purchase and maintenance of capital, fuel costs for transit between fishing grounds and ports as well as during gear operation, and labor. In the model, cost is simply proportional to effort:

$$\begin{array}{}\text{(9)}& {\mathrm{cost}}_k=c_k{E}_k\mathrm{d}t,\end{array}$$

where c_k is the cost per unit effort ($\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$).

When revenue exceeds costs, fishing effort in Eq. (6) increases. Any nonzero effort will lead to catches,

$$\begin{array}{}\text{(10)}& h_k\mathrm{d}t\mathrm{d}m=q_k{\mathit{\sigma }}_k{E}_k{f}_k\mathrm{d}t\mathrm{d}m,\end{array}$$

which couple the economic and biological modules in Eq. (1). A catch limit is imposed for numerical stability (i.e., to prevent harvesting more fish than the biomass available in each grid cell). The total catch rate for each group is then given by $C_k={\int }_{m_{\mathrm{0}}}^{m_{\mathrm{\infty },k}}{h}_k\mathrm{d}m$ (in $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). Note that when computing catches, but also costs and effort change, we set a lower limit on effort $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{50}}$ to allow the development of fishing and prevent division by zero in Eq. (6).

In BOATSv1, the ex-vessel fish price p_k is generally assumed to be constant in space and time ($\mathrm{1.264}\times {\mathrm{10}}^{-\mathrm{3}}$ USD g⁻¹), since observations suggest small historical variations (Sumaila et al., 2007; Galbraith et al., 2017). Similarly, cost c_k is also assumed to be constant ($\mathrm{1.852}\times {\mathrm{10}}^{-\mathrm{7}}$ $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$). Catchability q_k increases annually at a 5 % rate that accounts for sustained technological improvements and is the only temporally varying economic “forcing”. Empirical studies have estimated an average annual rate of 2 %–8 % between fisheries and periods. We select an annual rate of 5 % increase as, after testing when other observed economic parameters are forced, it accurately reproduces the historical development of fisheries with BOATS (Galbraith et al., 2017; Scherrer and Galbraith, 2020). A list of economic parameters and quantities is provided in Table 1; additional details are provided in Appendix A.

Table 1 List of parameters and quantities of BOATSv2. For a full list and values see Appendix A.

3 BOATSv2

Here we describe the features of BOATSv2 that provide an update to the original BOATSv1. Two of these features were introduced incrementally, in previously published work, in order to capture iron limitation in regions where iron is known to be scarce (Galbraith et al., 2019) and to represent management of fisheries (Scherrer and Galbraith, 2020; Scherrer et al., 2020). We provide a brief summary of these previous updates, before discussing the novel features added to the model in detail (see also Table 1 and Appendix A).

3.1 Previously published features

3.1.1 Reduced growth rates in iron-limited regions

Iron limitation of phytoplankton growth is widely recognized in the ocean, most prominently in HNLC regions (Tagliabue et al., 2017). Less is known about iron limitation of higher trophic levels in the ocean, including fish (Le Mézo and Galbraith, 2021). When satellite-based observational estimates of primary production are used as forcings, BOATSv1 overestimates fishing effort in HNLC regions, likely by simulating excessive biomass. Evidence of lack of adaptation by fish to low-iron regions suggests that low iron availability also significantly limits fish growth and could contribute to reducing fish abundance in large portions of the high seas (Galbraith et al., 2019).

Following Galbraith et al. (2019), we parameterize iron limitation of fish by reducing the trophic efficiency α, which determines the fraction of biomass incorporated into new tissues at each trophic step in HNLC regions:

$$\begin{array}{}\text{(11)}& {\mathit{\alpha }}^{\mathrm{v}\mathrm{2}}=\mathit{\alpha }\left({\displaystyle \frac{k_{{\mathrm{NO}}_{\mathrm{3}}^{-}}}{k_{{\mathrm{NO}}_{\mathrm{3}}^{-}}+{\mathrm{NO}}_{\mathrm{3}}^{-}}}\right).\end{array}$$

Here, the surface concentration of nitrate (NO${}_{\mathrm{3}}^{-}$, in µM) is taken as a proxy for iron limitation (Moore et al., 2013) and as an indicator of regions where fish are expected to be limited by the lack of iron, given the absence of other robust global estimates of surface iron concentrations or plankton iron contents. Note that here and in the following sections, the superscript “v2” indicates corrected quantities compared to the initial formulation in BOATSv1. This parameterization smoothly decreases the trophic efficiency as surface nitrate increases. The constant $k_{{\mathrm{NO}}_{\mathrm{3}}^{-}}=\mathrm{5}$ µM controls the strength of this effect and is constrained empirically. Nitrate concentrations are taken as the monthly minimum from the World Ocean Atlas (WOA) climatology (Locarnini et al., 2006) (see Appendix B). Although uncertainties persist regarding the impact of iron limitation on marine predators, this parameterization effectively reduces an overestimate in fishing effort in the Southern Ocean, North Pacific, equatorial Pacific, and, to some extent, North Atlantic in BOATSv1 (Galbraith et al., 2019).

3.1.2 Management with varying effectiveness

In BOATSv1, effort was generally assumed to follow an open-access dynamic equation (Eq. 6). This was modified to represent the influence of regulation by adjusting fishing effort to align with a prescribed target E_targ,k (Scherrer and Galbraith, 2020) as

$$\begin{array}{}\text{(12)}& {\left({\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{E}_k\right)}^{\mathrm{v}\mathrm{2}}=\left({\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{E}_k\right)e^{-S}+\left(\mathrm{1}-e^{-S}\right){\mathit{\kappa }}_s\left(E_{\mathrm{targ},k}-E_k\right),\end{array}$$

where S is a non-dimensional parameter representing how effectively the target is enforced. When S=0 the model follows open-access dynamics; when S>0 the nominal effort is nudged towards the target at a rate proportional to the regulation response parameter κ_s (${\mathrm{m}}^{-\mathrm{2}}{\mathrm{USD}}^{-\mathrm{1}}$).

This feature showed that strong fishery regulation is required to prevent overfishing if technological progress keeps increasing, making management effectiveness a key factor in future scenarios (Scherrer and Galbraith, 2020). For the rest of the paper, we set S=0 and focus on simulation of historical fisheries up to the time they reached maximum catches, for which the open-access dynamic was shown to be adequate (Galbraith et al., 2017; Guiet et al., 2020).

3.2 Newly added features

3.2.1 Separate pelagic and demersal energy pathways

Variations in energy input at the base of marine food webs significantly affect biomass accumulation and cycling, thereby altering the sensitivity of different fish communities to climate and environmental factors (Petrik et al., 2019). Pelagic communities are more tightly tied to surface planktonic production (i.e., net primary production Π_ψ), whereas benthic communities depend on the delivery of organic material to the seafloor (i.e., particle flux at bottom Π_β in $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) (Stock et al., 2017; van Denderen et al., 2018). The two types of communities also experience different temperatures, with surface temperature (here, T=T₇₅) controlling the metabolic rates of pelagic fish and bottom temperature (T=T_bot) that of demersal fish.

To account for these ecological differences, we modified BOATS to resolve separate pelagic and demersal fish communities. Both communities are described by the same set of equations described above, Eqs. (1)–(12), but solved separately with independent sets of environmental forcings (see environmental forcing in Fig. 1). Pelagic fish are forced by surface conditions (Π_ψ and T=T₇₅), while demersal fish are forced by bottom conditions (Π_β and T=T_bot). Whereas the energy supply to the pelagic community remains dependent on surface net primary production (NPP), the particle flux to the bottom provides the energy input to the demersal community. The particle flux is modeled as a depth-dependent fraction of surface primary production:

$$\begin{array}{}\text{(13)}& {\mathrm{\Pi }}_{\mathit{\beta }}={\mathrm{\Pi }}_{\mathit{\psi }}\times {\mathrm{pe}}_{\mathrm{ratio}}\times {\left({\displaystyle \frac{z_{\mathrm{bot}}}{z_{\mathrm{eu}}}}\right)}^{b_{\mathrm{a}}}.\end{array}$$

This formulation assumes a power-law decrease in the flux of organic material below the euphotic layer, i.e., a typical Martin curve (Martin et al., 1987; Buesseler and Boyd, 2009). The attenuation coefficient $b_{\mathrm{a}}=-\mathrm{0.8}$ is selected within the range of plausible values (Gloege et al., 2017), and the euphotic layer depth z_eu=75 m is assumed to be fixed, although both could be modeled to vary in space and time. Here, similar to prior work (Stock et al., 2017; van Denderen et al., 2018; Petrik et al., 2019), we focus on first-order variations in fish biomass in increasingly deep habitats, where food becomes progressively scarce. We calculate the term $(z_{\mathrm{bot}}/z_{\mathrm{eu}}{)}^{b_{\mathrm{a}}}$ using the high-resolution bathymetry $z_{\mathrm{bot}}^{\prime }$ from the ETOPO global surface relief at $\mathrm{1}/\mathrm{10}$th degree (Amante and Eakins, 2009), taking the average across 1° grid cells $(z_{\mathrm{bot}}/z_{\mathrm{eu}}{)}^{b_{\mathrm{a}}}=\stackrel{\mathrm{‾}}{(z_{\mathrm{bot}}^{\prime }/z_{\mathrm{eu}}{)}^{b_{\mathrm{a}}}}$. Note that when $z_{\mathrm{bot}}^{\prime }$ is shallower than z_eu, export is only determined by pe_ratio, which is taken as a function of surface temperature T₇₅ (in °C) and net primary production Π_ψ, following Dunne et al. (2005).

In the pelagic ocean, the typical size of phytoplankton, m_ψ, varies markedly between productive and oligotrophic regions. This variation affects both the length of the food web and the proportion of production accessible to fish communities (Ryther, 1969). We use an empirical phytoplankton size model to account for this variation (Dunne et al., 2005). In analogy with the pelagic ecosystem, we assume that the representative size of benthic organisms at the base of the demersal food web, m_β, influences the fraction of energy that reaches demersal fish. For simplicity, we take m_β to be globally uniform. Unlike m_ψ, for which empirical parameterizations exist, m_β is poorly constrained. We keep most food web parameters the same for pelagic and demersal fish, with the exception of the activation energy for growth ω_a,A and mortality ω_a,λ, since observations of growth rates suggest significant differences between the two communities (van Denderen et al., 2020).

3.2.2 Heterogeneous costs

Simulations with BOATSv1 suggest that variations in the cost per unit effort of fishing c_k (in $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) played only a secondary role in the development of global fisheries (Galbraith et al., 2017). Yet, heterogeneous costs in the global ocean could modulate the spatial distribution of fishing effort and its evolution over time (Swartz et al., 2010; Anticamara et al., 2011; Lam et al., 2011). Reconstructions of fishing effort in the high seas suggest more than 2-fold average cost differences between distinct fishing gears and regions (see Sala et al., 2018, Kroodsma et al., 2018, and Appendix C).

To simulate the effect of heterogeneous fishing costs on the historical offshore expansion of fisheries, we replaced the constant costs per unit effort (c_k=5.85 $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) in BOATSv1 by spatially varying costs using a linear function of the distance to shore for effort targeting pelagic fish, d_coast (in km), and a linear function of seafloor depth z_bot (in m) for effort targeting demersal fish (see economic forcing Fig. 1):

$$\begin{array}{}\text{(14)}& c_k^{\mathrm{v}\mathrm{2}}\left(x=d_{\mathrm{coast}},z_{\mathrm{bot}}\right)=\left\{\begin{array}{l}{c}_k\mathrm{when}x\le x^{*}\\ c_k+\mathit{\delta }(x-x^{*})\mathrm{when}x>x^{*},\end{array}\right.\end{array}$$

where x^* is a reference parameter that determines the boundary between coastal and high seas regions. For pelagic effort, x^* identifies a coastal band over which transit costs are assumed to be small compared to other costs. Here we adopt $x^{*}=\mathrm{370}$ km (or 200 nm), the limit of exclusive economic zones separating coastal regions and high seas. For demersal effort, x^* identifies a depth threshold above which the cost of setting and hauling gears is negligible compared to other costs, and we set $x^{*}=\mathrm{200}$ m. The parameter δ is the proportionality constant for the increase in costs beyond these coastal bands (in $\mathrm{USD}{\mathrm{km}}^{-\mathrm{1}}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for pelagic effort and $\mathrm{USD}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for demersal effort).

For distance-dependent costs, we select $\mathit{\delta }=\mathrm{7.9}\times {\mathrm{10}}^{-\mathrm{3}}$ $\mathrm{USD}{\mathrm{km}}^{-\mathrm{1}}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ such that the average high seas fishing cost is 9.3 $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, comparable with an empirical upper mean value of 8.9 $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (see Appendix C). For depth-dependent costs, the depth of the fishing grounds z_bot is determined from high-resolution ($\mathrm{1}/\mathrm{10}$th degree) bottom topography (Amante and Eakins, 2009), taking the shallowest depth in each 1° model grid cell. We set $\mathit{\delta }=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$ $\mathrm{USD}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ such that the average high seas fishing cost is 9.9 $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for exploitation of deep demersal stocks, comparable with an empirical upper boundary of mean high seas trawling costs of 9.2 $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Appendix C). Given the uncertainty regarding whether distance or seafloor depth has a greater impact on costs in pelagic and demersal fisheries, we first tested distance- and depth-dependent costs separately and then added them to determine their combined impact.

3.2.3 Heterogeneous catchability

In BOATSv1, technological progress, represented by an exponential increase in the catchability q_k at a rate of 5 % yr⁻¹, was shown to be a dominant driver of the development of fisheries (Galbraith et al., 2017). While a homogeneous increase rate approximates the first-order effect of technological progress well, heterogeneous technological efficiencies among fisheries could modulate this development across regions (Palomares and Pauly, 2019), especially as separate gears target distinct resources and are deployed in different ecosystems (Kroodsma et al., 2018). Similar to cost, spatially heterogeneous catchability could have influenced the spatial expansion of fisheries and/or the deepening of catches with time (Watson and Morato, 2013).

To simulate the effect of heterogeneous catchability, the exponential increase is spatially weighted (see economic forcing in Fig. 1):

$$\begin{array}{}\text{(15)}& q_k^{\mathrm{v}\mathrm{2}}\left(x=z_{\mathrm{bot}},y=\mathrm{P}\mathrm{or}\mathrm{D}\right)=q_k\mathrm{Pr}\left(x\right)\mathrm{Of}\left(y\right),\end{array}$$

where Pr(x=z_bot) accounts for spatial variations in technological efficiencies with seafloor depth, and Of(y=P or D) is an offset between the catchability of pelagic and demersal resources.

Commercially exploited fish often aggregate near seamounts and other shallow features, resulting in the local establishment of fisheries (e.g., 57 % of longlining activity; Kerry et al., 2022). The coarse resolution of BOATS prevents a direct representation of seamounts. However, the presence of seamounts could increase both fish biomass density and profitability within a model's grid cell, as opposed to the case where resources were more homogeneously distributed across the grid cell. The profile Pr(x=z_bot) parameterizes the effect of seamounts and more generally an increase in the density of resources in shallow regions:

$$\begin{array}{}\text{(16)}& \mathrm{Pr}(x=z_{\mathrm{bot}})=q_{\min }+(\mathrm{1}-q_{\min }){\displaystyle \frac{{\log }_{\mathrm{10}}\left(x_{\max }\right)-{\log }_{\mathrm{10}}\left(x\right)}{lo{g}_{\mathrm{10}}\left(x_{\max }\right)-{\log }_{\mathrm{10}}\left(x_{\mathrm{mean}}\right)}}.\end{array}$$

Here, q_min=0.8 is the minimum efficiency of gears targeting pelagic resources (see Appendix D), and x_mean and x_max at 2372 and 5750 m depth, respectively, are based on the median and deepest depths of seamounts where fishing activity occurs (Kerry et al., 2022). The depth of the fishing grounds z_bot is determined from ETOPO at $\mathrm{1}/\mathrm{10}$th resolution (Amante and Eakins, 2009), coarsened by taking the shallowest depth in each 1° model grid cell, as described above.

The dominant gears used to target different communities (pelagic vs. demersal) are characterized by different efficiencies (see Appendix D). We tested the effect of separate catchabilities for pelagic and demersal communities, setting Of(P)=1.4 when Of(D)=1. This offset was estimated from the technology coefficients of different gear targeting pelagic and demersal species, weighted by the fraction of global fishing effort for 16 different gears (Kroodsma et al., 2018; Palomares and Pauly, 2019).

3.3 Forcing and initialization

Forcing BOATSv2 requires surface temperature T₇₅, bottom temperature T_bot, and net primary production Π_ψ. Since we are interested in the recent ocean state, we use climatological observations, and, to assess improvements between BOATS versions, we adopt the same forcing as in Carozza et al. (2017). Surface temperature (T₇₅) and temperature at the seafloor (T_bot) are taken from the World Ocean Atlas 2009 (Locarnini et al., 2006). T₇₅ (in °C) is calculated as the mean temperature over the top 75 m on a 1° grid. T_bot (in °C) is calculated by averaging temperatures at different depths, weighted by the fraction of each depth within a model grid cell as reported by the ETOPO $\mathrm{1}/\mathrm{10}$ bathymetry dataset (Amante and Eakins, 2009). Recognizing that the resolution of observational temperature datasets such as the WOA decreases with depth, we select the layers closest to the bottom as indicative of the temperature near the seafloor. For Π_ψ we take the average of three satellite-based estimates at 1° resolution (Behrenfeld and Falkowski, 1997; Carr et al., 2006; Marra et al., 2007). Note that Π_ψ and then Π_β are forced once converted to $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. Forcing BOATS with two-dimensional grids does not account for vertical positions along the water column but characterizes mean environmental conditions where many harvested fish live.

The model is initialized by a “pristine” ocean biomass distribution determined by spinning up the model without fishing for 300 years to reach a steady state. Then, the ecological and economic modules are run together with an increase in the catchability q_k for another 300 years, starting with a small, globally uniform effort. In a given region, fishing begins once catches become profitable (i.e., revenue_k>cost_k in Eq. 6). The open-access dynamics generally first drive an increase in catch, followed by a peak and decline due to overfishing (Guiet et al., 2020). To align simulations with observations (see Sect. 4.1), we estimate the time of the peak catch integrated across LMEs and align it with the time of the observed peak catch, which occurs in the 1990s (Pauly and Zeller, 2016).

4 Parameterization procedure

In this section, we first describe the observations used for the evaluation of BOATSv2 and then detail the procedure used to parameterize the model, which is based on the following two steps (Fig. 2). (1) Ecological update: we start by focusing on coastal regions, where most of the observed catch originates and where economic parameterizations are more homogeneous. We parameterize separate pelagic and demersal pathways (Π_ψ and Π_β) and growth limitation in HNLC regions for pelagic species (α^v2), to determine the best parameter sets for 11 undetermined parameters of the ecological module (see Table 2). (2) Economic update: we then fine-tune the parameters of the economic module using spatially heterogeneous costs and catchability, considering the global ocean (i.e., $c_k^{\mathrm{v}\mathrm{2}}$ and $q_k^{\mathrm{v}\mathrm{2}}$). Results of this parameterization procedure are described in Sect. 5.1.

Figure 2 Schematic diagram of the parameterization procedure starting from BOATSv1 (in blue), with two steps. (1) Ecological update (in green): a Monte Carlo tuning procedure with five selection criteria is applied to a modified version of BOATSv1 that allows separate pelagic and demersal pathways and growth limitation in iron-limited regions. An ensemble of 20 000 simulations is carried out only for coastal regions with various parameter sets, and we identify a set of 5 (or 10 extended) best parameter sets. (2) Economic update (in red): with three selection criteria, we identify the best economic parameterizations applied to the optimized intermediate BOATSv2-Bio version to determine BOATSv2. We use simulations that include the high seas for the five best parameter sets. Observations used for the parameterization of both are shown in gray.

4.1 Observational data and diagnostics for model evaluation

We use multiple empirical data sources, including catch, biomass, and fishing effort, to tune and evaluate BOATSv2 (in gray Fig. 2). Comparisons are made on globally integrated quantities, quantities integrated across LMEs to assess regional variability in coastal regions, and quantities integrated beyond the boundary of LMEs, i.e., across high seas ecosystems (HSEs; see Appendix E), to assess variability in the open ocean away from coastal influences. We focus on observations around the peak catch in the 1990s but also include observations in the 1950s and 2000s for additional insight.

For fish catch, we use two catch reconstruction datasets: (1) the first is from The Sea Around Us project (SAU; Pauly et al., 2020), corrected for under-reported catch. For the SAU, catches by functional type allow separating pelagic (P) and demersal (D) species (see Appendix F). (2) The second is the database from Watson (2017) (hereafter WAT), including wild catch and corrected for illegal and unreported fisheries. When comparing catches by LME, we focus on 55 LMEs (out of 66) and ignore the Black Sea and a number of high-latitude regions to avoid errors caused by biases in satellite-based chlorophyll and the lack of representation of the effects of sea ice on the marine ecosystem (as in Carozza et al., 2017). We define two diagnostics to help correct for biases in BOATSv1: the fraction of catch in the high seas, $R_{\mathrm{C}}=C_{\mathrm{HS}}/(C_{\mathrm{HS}}+C_{\mathrm{CS}})$, and the catch-weighted mean depth of fishing $Z_{\mathrm{C}}=\left({\sum }_{\mathrm{lat},\mathrm{lon}}C{z}_{\mathrm{bot}}\right)/\left({\sum }_{\mathrm{lat},\mathrm{lon}}C\right)$ (in m).

For biomass observations, we use the RAM Legacy Stock Assessment Database (Ricard et al., 2012). Stock assessment data are used to estimate mean catch to biomass ratios (C:B hereafter) in 25 LMEs where enough stock assessments are available, following Bianchi et al. (2021). We also compare historical changes in fish biomass to a global reconstruction based on stock assessments (Worm and Branch, 2012). Furthermore, we compare the model with two biomass databases derived from fishery-independent surveys: the first, encompassing demersal species across 14 large marine ecosystems (LMEs) in the Northern Hemisphere, ranging from the Bering Sea to northern Europe, is based on a recent synthesis of bottom trawl data (van Denderen et al., 2023; Maureaud et al., 2024). The second, focusing on pelagic species, is built on standardized trawls of coastal pelagic species in two LMEs along the North American west coast (Zwolinski et al., 2012).

Finally, we include a comparison with reconstructions of nominal effort for the global fishing fleet in both artisanal and industrial sectors to shed light on the regional development of fisheries (Rousseau et al., 2019, 2024). Similar to catch, we focus on a subset of 55 LMEs where model forcings are more suitable.

4.2 Ecological update: Monte Carlo ensemble

To calibrate the revised model based on BOATSv1 (in blue Fig. 2), we must specify the value of 11 poorly determined parameters (Table 2). (Note that the trophic scaling is a function of two free parameters and is thus completely determined by their values.) These parameters are not well constrained by the literature: the activation energies ${\mathit{\omega }}_{\mathrm{a},\mathrm{A}-\mathit{\lambda }}$; the scaling exponents b and h; the constants A₀ and ζ₁; the trophic scaling $\mathit{\tau }=\ln \left(\mathit{\alpha }\right)/\ln \left(\mathit{\beta }\right)$ (Brown et al., 2004), itself a function of trophic efficiency α and predator–prey mass ratio β (Barnes et al., 2010); the egg survival fraction s_e; the threshold mass for fishing selectivity $e_{m_{\mathrm{\Theta },k}}$; and the representative size of benthic organisms m_β. Following previous work (Carozza et al., 2017; Bianchi et al., 2021), we adopt an ensemble Monte Carlo approach, running replicates with parameter sets randomly chosen from plausible ranges of values (see Table 2). Note that compared to the calibration of BOATSv1 (Carozza et al., 2017), we assume a subset of parameters (k_E, Π^*, and c_σ) to be relatively well constrained, since previous analysis showed that variation in these parameters had no significant effect. We also updated a few prior parameter ranges based on recent analyses like α (Stock et al., 2017; Eddy et al., 2020) and because previously optimized values were close to the boundaries of the ranges (b, s_e). We keep the same parameters for both pelagic and demersal communities, except for the temperature dependence of growth and mortality.

We run 20 000 simulations, each with a distinct combination of parameters, integrated with gradually increasing catchability over time, and select the best simulations according to global and local criteria. These criteria are updated from Carozza et al. (2017) to provide an evaluation of the model performance in reproducing the following features of pelagic and demersal communities.

1. Global catch. The best simulations must predict observed global fish catch, when integrated over the 55 LMEs, for the years of maximum catch in the 1990s, i.e., $C_{\mathrm{Globmax}}^{\text{SAU-WAT}}\simeq \mathrm{100}\times {\mathrm{10}}^{\mathrm{6}}$ t yr⁻¹ (using tons of wet biomass). There are significant uncertainties around these reconstructions; furthermore, migratory species not represented in BOATS can influence the model's maximum yields. Therefore, we allow catches to be within the range $C_{\mathrm{Globmax}}\in [\mathrm{70},\mathrm{150}]\times {\mathrm{10}}^{\mathrm{6}}$ t yr⁻¹.

2. Demersal to pelagic catch. At the global catch peak, simulations must capture the fraction of pelagic and demersal catch. Integrated over all 55 LMEs, pelagic catch and demersal catch respectively account for 45 % and 55 % of catches in SAU – that is, the ratio between demersal and pelagic catch at the peak is $R_{\mathrm{Globmax},\mathrm{D}/\mathrm{P}}^{\mathrm{SAU}}\simeq \mathrm{1.2}$. Because of uncertainties around the SAU reconstructions, the presence of migratory species, and additional uncertainty in the allocation of pelagic vs. demersal catch (Appendix F) we allow this ratio to vary within the range $R_{\mathrm{Globmax},\mathrm{D}/\mathrm{P}}\in [\mathrm{0.8},\mathrm{1.8}]$.

3. Catch to biomass. To ensure that global catches are supported by realistic rates of fish biomass production, we compare the model catch to biomass ratio (C:B) averaged over 25 LMEs to the observational estimate from the RAM database (see Bianchi et al., 2021). We retain simulations for which a Kolmogorov–Smirnov test indicates that the modeled C:B ratios follow the same distribution as the stock assessment data, rejecting cases where distributions are found to be different at the 1 % significance level.

4. Size abundance. To preserve a realistic partitioning of fish catches by size groups, for the best simulations we constrain the catch of medium and large sizes to be in the observed range relative to fish in the small group, i.e., 0.3 C_small<C_med and $\mathrm{0.1}{C}_{\mathrm{small}}<C_{\mathrm{lrg}}<\mathrm{0.8}{C}_{\mathrm{small}}$.

5. Spatial variability near the coast. Finally, we assess the regional variability of catch at the time of the global peak by computing Pearson correlation coefficients of simulated catch densities compared to observations ($r_{\mathrm{LME}\mathrm{90}\mathrm{s}}^{\mathrm{SAU}}$ or $r_{\mathrm{LME}\mathrm{90}\mathrm{s}}^{\mathrm{WAT}}$). We also compare simulated and observed maximum catch per functional type and LME, independently of the peak year, to estimate the model capability to reproduce maximum yields per group ($r_{\mathrm{LMEmax}}^{\mathrm{SAU}\mathrm{P}}$ or $r_{\mathrm{LMEmax}}^{\mathrm{SAU}\mathrm{D}}$).

Criteria (1) to (5) identify parameter sets (Table 2 and results Sect. 5.1.1) that best capture global properties of catches and of fish production per unit biomass for both pelagic and demersal species, focusing on well-sampled coastal regions, completing the ecological update to BOATSv2-Bio (in green in Fig. 2).

Table 2 Model parameters and summary results for the Monte Carlo ensemble (update of Table 3 in Carozza et al., 2017). The sampling distribution of each parameter used in the Monte Carlo simulation are shown, where X(p1,p2) represents the probability distribution (N for normal, U for uniform), and p1 and p2 are the mean and standard deviation of each parameter, respectively, for the pelagic P and demersal D communities when it applies. Opt. refers to the subset of optimized Monte Carlo simulations, N.O. to the subset of non-optimized simulations. SD is the standard deviation, and KS p value is the p value of the two-sample Kolmogorov–Smirnov test applied to the optimized and non-optimized sets. The three last variables are fixed compared to previous optimizations.

^* We merge two ensembles of 10 000 simulations each, with slightly different distributions for b. The first ensemble prompted re-selection of the parameter range for the second.

4.3 Economic update: sensitivity to cost and catchability

After improving the ecological realism of the simulations by tuning selected parameters with a Monte Carlo approach, we improve the economic realism by incorporating heterogeneous costs ($c_k^{\mathrm{v}\mathrm{2}}$) and catchability ($q_k^{\mathrm{v}\mathrm{2}}$) (in red in Fig. 2). We further evaluated the effect of considering different economic parameterizations and selected the best combination based on regional and global criteria. This evaluation compares the following.

• 6.

  High seas variability. Beyond LMEs, in HSEs, once heterogeneous costs or catchability are activated we compute Pearson correlation coefficients at peak r_HSE90s to weigh improvements for predicted catch, similarly to step (5) in coastal regions.

• 7.

  High seas catch fraction. This constraint determines how costs and catchability influence the catches in the high seas and indirectly the historical offshore expansion of fisheries. We computed the high seas catch fraction in the 1950s (${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{50}\mathrm{s}}$) and near the global catch peak of the 1990s (${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{90}\mathrm{s}}$). We expect ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{50}\mathrm{s}}\simeq \mathrm{0.06}$ to increase to only ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{90}\mathrm{s}}\simeq \mathrm{0.09}$ at the global peak, while catch fractions for pelagic and demersal fish increase from ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{50}\mathrm{s}}=\mathrm{0.10}$ to ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{90}\mathrm{s}}=\mathrm{0.11}$ and ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{50}\mathrm{s}}=\mathrm{0.05}$ to ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{90}\mathrm{s}}=\mathrm{0.07}$, respectively.

• 8.

  Mean fishing depth. Finally, to better characterize the offshore expansion of fisheries, we computed the catch-weighted mean depth over which fishing occurs in the 1950s and 1990s, ${\stackrel{\mathrm{‾}}{Z}}_{\text{C50s–90s}}$. For demersal catch (${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{50}\mathrm{s},\mathrm{D}}<\mathrm{136}$ m and ${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s},\mathrm{D}}<\mathrm{206}$ m), the mean fishing depth reflects the historical deepening of fishing grounds (Watson and Morato, 2013). For pelagic catch (${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{50}\mathrm{s},\mathrm{P}}<\mathrm{266}$ m and ${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s},\mathrm{P}}<\mathrm{546}$ m), it reflects an offshore expansion of fishing effort towards high seas regions with a deeper seafloor.

Criteria (6) to (8) allow refining our understanding of the regional variability of catch and their sensitivity to economic parameterizations (see results in Sect. 5.1.2). They are applied to simulation sets with optimum parameter sets, for multiple cost and catchability profiles, to best capture regional properties of catches, ultimately determining BOATSv2 (Fig. 2).

5 Results and discussion

5.1 Parameterization

5.1.1 Ecological parameters

The Monte Carlo ensemble allows the identification of optimum ecological parameter sets (Fig. 2). All simulations span a total catch range of more than 6 orders of magnitude (see Fig. 3). We find that 12 % of simulations satisfy the first criterion of global catch (see Sect. 4.2). Then, for criterion (2), the ratio of demersal to pelagic catch varies by more than 3 orders of magnitude, with 20 % of the simulations with realistic ratios, leaving us with 3.0 % when combined with criterion (1) (Fig. 3). Among these, simulations capturing the observed catch to biomass ratio leave us with 0.8 % of all simulations, and acceptable size distributions ultimately lead to 0.2 % (42 simulations) of all simulations satisfying criteria (1) to (4).

Figure 4 shows the time series of catch, nominal effort, and biomass over all 55 LMEs from 1900 to 2050 for the 42 simulations that meet criteria (1) to (4). In each simulation, global catch increases until reaching a peak, beyond which biomass depletion limits recruitment and drives a fall in catch. Effort follows a comparable pattern, but with a consistent time lag. Biomass monotonically decreases from an initial, nearly pristine state. These features are comparable to observational reconstructions (Fig. 4a, C^SAU and C^WAT). A delayed response of nominal effort is also consistent with observations (Fig. 4b), while the consistent decrease in biomass compares well with aggregated stock assessment data normalized to the pristine period (Fig. 4c, d).

Similar to prior work with BOATSv1, we focus the next analysis on the five best ensemble members selected to capture parameter uncertainty, while maintaining reasonable computational costs. These five parameter sets are selected out of the 10 best of the 42 simulations based on their ability to reproduce regional variability in peak catches by LMEs – criterion (5). The peak catch is determined almost exclusively by ecological parameters, making it a valuable way to discriminate amongst them (Carozza et al., 2017). Accordingly, we rank the 42 simulations by the Pearson correlation coefficient of simulated vs. observed catch in the 55 LMEs ($r_{\mathrm{LMEmax}}^{\mathrm{SAU}}$; see Fig. 3) and select 5 ensemble members out of the top 10. These five chosen parameter sets comprise diverse shapes of catch, effort, and biomass histories, but, once averaged together, they provide an ensemble mean that matches the observed historical development of these quantities across LMEs ($r_{\mathrm{LMEmax}}^{\mathrm{SAU}}\in [\mathrm{0.63},\mathrm{0.69}]$; see Table 3 and Fig. 4, dark red lines).

Note that the Pearson correlation coefficients r_LME90s between observed and simulated catch by LME at global peak are comparable with and without updated ecological features (see BOATSv1 compared to other model variants; Table 4; $r_{\mathrm{LME}\mathrm{90}\mathrm{s}}^{\mathrm{SAU}}\simeq \mathrm{0.69}$ and $r_{\mathrm{LME}\mathrm{90}\mathrm{s}}^{\mathrm{WAT}}\simeq \mathrm{0.73}$). However, the updated ecological features provide large improvements in the high seas (e.g., $r_{\mathrm{HSE}\mathrm{90}\mathrm{s}}^{\mathrm{SAU}}$ increases from 0.22 in v1 to 0.58). This improvement in the high seas is partly explained by the representation of iron limitation on fish growth ($r_{\mathrm{HSE}\mathrm{90}\mathrm{s}}^{\mathrm{SAU}}$ increases to 0.81 from 0.22 in v1), while along coastal regions, iron limitation alone is insufficient to explain catch (Table 4).

Figure 3 Simulations of global catch features from the BOATSv2 Monte Carlo ensemble. The total catch at the global peak of the 1990s C_Globmax is shown as a function of the ratio between demersal and pelagic catch at the global catch peak, $R_{\mathrm{Globmax},\mathrm{D}/\mathrm{P}}$, for the 20 000 simulations in the ensemble. Colors show the Pearson correlation coefficient $r_{\mathrm{LMEmax}}^{\mathrm{SAU}}$ of simulated vs. observed (SAU) maximum catch in 55 LMEs, C_LMEmax, for each simulation. The dotted black horizontal line shows the reference global harvest, $C_{\mathrm{Globmax}}^{\text{SAU-WAT}}=\mathrm{100}\times {\mathrm{10}}^{\mathrm{12}}$ g yr⁻¹. The dotted black vertical line shows the observed ratio of demersal to pelagic catch at the global peak, $R_{\mathrm{Globmax},\mathrm{D}/\mathrm{P}}^{\mathrm{SAU}}=\mathrm{1.2}$. The horizontal and vertical gray lines indicate the ranges within which the best simulations are selected.

Figure 4 Historical simulations of catch, effort, and biomass for the best 42 simulations in the BOATSv2 Monte Carlo ensemble, forced only with exponentially increasing catchability over time. (a) Catch, (b) nominal effort, (c) biomass, and (d) biomass normalized to initial biomass in the selected 55 LMEs for the 42 best parameter sets from 1900 through 2050. Ensembles are aligned at the catch peak of the 1990s. The light blue lines show each parameter set; the dark blue lines show the 10 best simulations, out of which the 5 purple lines show the final best ensemble. The thick red line is the mean of the five best ensemble simulations. Black lines in panels (a), (b), and (d) show observational reconstructions, consisting of catch from SAU (black line in panel a, C^SAU) and the WAT database (dotted black line in panel a, C^WAT), effort (black line in panel b), and biomass from fish stock assessments normalized to the initial state (black line in panel d). Note that the simulations do not include the effects of climate change, environmental variability other than the seasonal cycle, and management.

Table 3 List of parameter values for the selected five best simulations with updated ecological features.

Table 4 Sensitivity of model variants for various combinations of ecological and economical processes. We report Pearson coefficients between mean simulated and observed catch across the 55 selected LMEs at the global peak r_LME90s or for the historical maximum r_LMEmax and across high seas ecosystems (HSEs) at global the peak r_HSE90s. We report the fraction of high seas vs. coastal catch R_C and the mean depth of catch ${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}}$ in the 1950s and near the global peak of the 1990s. The superscripts SAU and WAT indicate the kind of observation for pelagic (P) or demersal (D) communities with SAU.

^* $\text{v2-Bio}=v\mathrm{1}+{\mathit{\alpha }}^{\mathrm{v}\mathrm{2}}+{\mathrm{\Pi }}_{\mathit{\psi }}$.

5.1.2 Economic parameters

As summarized in Table 4, increased ecological realism improves the model's ability to reproduce high seas fisheries, in particular the fraction of high seas catch (${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{90}\mathrm{s}}$ down to 0.16 from 0.40 in BOATSv1) and the catch-weighted mean depth of fishing (${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s}}$ down to 694 from 1698 m). These improvements reflect growth limitation in HNLC regions (α^v2) (Galbraith et al., 2019) and, to a greater extent, explicit separation of pelagic and demersal energy pathways (Π_ψ and Π_β+α^v2; hereafter BOATSv2-Bio; green line in Fig. 5a, b).

Since economic drivers could explain additional spatial variability, we test plausible heterogeneous economic parameterizations (Fig. 2). The heterogeneous costs and catchabilities have little impact on the coastal variability of catch at the time of the global peak ($r_{\mathrm{LME}\mathrm{90}\mathrm{s}}^{\mathrm{SAU}}\simeq \mathrm{0.69}$ and $r_{\mathrm{LME}\mathrm{90}\mathrm{s}}^{\mathrm{WAT}}\simeq \mathrm{0.73}$; Table 4). Note that the comparison reveals better correlations when comparing models with WAT catch reconstructions instead of SAU reconstructions. Most of the improvement is explained by higher mean catches in Australian LMEs (compare Fig. 5c and d), but the explanation for such discrepancy in the observational reconstructions remains unclear. We also compare the Pearson correlation coefficients of maximum catch by LME for pelagic and demersal catch separately. Heterogeneous costs or catchability show no effect on the variability of maximum pelagic and demersal catch yields: $r_{\mathrm{LMEmax}}^{\mathrm{SAUP}}\simeq \mathrm{0.46}$ vs. $r_{\mathrm{LMEmax}}^{\mathrm{SAUD}}\simeq \mathrm{0.69}$ (Table 4); these should instead influence the timing of the development of fisheries. Both correlations suggest that, along the coast, catches are independent of economic parameterizations and are instead controlled mainly by the environment.

To select economic parameterizations, criterion (6) and the catch variability in the high seas indicate only minor variations; economic parameters do not significantly enhance the accuracy of catches in the HSEs compared to the improvement from BOATSv1 to BOATSv2-Bio (see Table 4). However, for criterion (7), economic parameterizations further reduce the high seas catch fraction (down to ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{90}\mathrm{s}}=\mathrm{0.10}$ compared to ≃0.09; Table 4), and criterion (8) significantly reduces the mean fishing depth (down to ${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s}}=\mathrm{315}$ m; Table 4).

Heterogeneous costs and catchability parameterizations have unequal effects. Depth- and distance-dependent costs respectively reduce the offshore expansion of demersal (${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s},\mathrm{D}}$ down to 148 m) and pelagic (${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s},\mathrm{P}}$ down to 873 m) catches. Since both slightly improve aspects of the simulations, we retain both parameterizations in the final BOATSv2 update. Regarding catchability, accounting for the effect of seamounts also reduces the development of fishing over a deep seafloor. When applied to pelagic catches only, ${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s}}$ decreases to 526 m (Table 4). When this correction is applied both to pelagic and demersal communities, this is further reduced to ${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s}}=\mathrm{315}$ m, close to observational estimates. However, the maximum depth of the demersal catch becomes excessively shallow (${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s},\mathrm{D}}=\mathrm{103}$ m). Therefore, we retain the heterogeneous catchability parameterization only for pelagic fishing effort in the final BOATSv2 update.

For the final reference simulations with BOATSv2, the inclusion of costs that increase with the distance from shore, costs that increase with seafloor depth for demersal fishing, and catchabilities that decrease with seafloor depth for the pelagic community leads to an improvement of the model's ability to reproduce the delayed development of high seas fisheries and the progressive deepening of catch (Table 4 and see red lines in Fig. 5a, b).

In summary, at the LME level, BOATSv2 and BOATSv1 have similar accuracy in their representation of regional catches (see r_LME90s=0.69 or 0.73 when compared to SAU or WAT; Table 4 and Fig. 5c, d). This lack of improvement in the new model version is explained by limited accuracy in predicting pelagic catches across LMEs (see $r_{\mathrm{LMEmax},\mathrm{P}}=\mathrm{0.46}$ vs. $r_{\mathrm{LMEmax},\mathrm{D}}=\mathrm{0.69}$; Table 4 and Fig. 5e, f). Nevertheless, BOATSv2 better captures the large-scale variability of catches in the HSEs, which are approximately 1 order of magnitude smaller than in LMEs (see r_HSE90s=0.51 or 0.64 compared to SAU or WAT; Table 4 and Fig. 5g, h), and better reproduces their historical offshore expansion (${\stackrel{\mathrm{‾}}{R}}_{\mathrm{C}\mathrm{90}\mathrm{s}}=\mathrm{0.11}$ and ${\stackrel{\mathrm{‾}}{Z}}_{\mathrm{C}\mathrm{90}\mathrm{s}}=\mathrm{420}$ m; Table 4).

Figure 5 Evaluation of BOATSv2. (a) Observed and simulated historical development of high seas vs. coastal seas catch R_C. (b) Observed and simulated historical deepening of global catches Z_C. (c, d) Scatter plot of the observed vs. simulated catch at the global peak in 55 LMEs for SAU reconstructions $C_{\mathrm{LME}\mathrm{90}\mathrm{s}}^{\mathrm{SAU}}$ in (c) and WAT reconstructions $C_{\mathrm{LME}\mathrm{90}\mathrm{s}}^{\mathrm{WAT}}$ in (d). (e, f) Scatter plot of the observed vs. simulated maximum catch in LMEs for pelagic $C_{\mathrm{LMEmax}}^{\mathrm{SAUP}}$ in (e) and demersal catch $C_{\mathrm{LMEmax}}^{\mathrm{SAUD}}$ in (f). (g, h) Scatter plot of the observed vs. simulated catch at the global peak in HSEs for SAU reconstructions $C_{\mathrm{HSE}\mathrm{90}\mathrm{s}}^{\mathrm{SAU}}$ in (g) and WAT reconstructions $C_{\mathrm{HSE}\mathrm{90}\mathrm{s}}^{\mathrm{WAT}}$ in (h). In panels (a) and (b), SAU and WAT reconstructions are indicated by solid and dotted lines, respectively; BOATSv1 is indicated by the blue line, BOATSv2-Bio by green, and BOATSv2 by red. In panels (c)–(h) blue dots and lines show BOATSv1 and red dots and lines BOATSv2. Numbers next to each dot indicate the LMEs or HSEs (see regions in Appendix E).

5.1.3 Model sensitivity

The best parameter sets selected by the Monte Carlo approach (Sect. 4.2) provide insights on the functioning of ecological communities. Of the 11 parameters that were optimized for, 6 have posterior distributions significantly different from the prior distributions (p values$<{\mathrm{10}}^{-\mathrm{2}}$; Table 2). The posterior distributions for these six parameters were also different when optimizing BOATSv1 (see Carozza et al., 2017), confirming their essential role in influencing the sensitivity of the model.

First, the ensemble mean allometric scaling exponent of b=0.72 (the range is 0.70–0.80 for the five best ensembles; see Table 3) is larger than the BOATSv1 value of 0.65 (Carozza et al., 2017) but in the middle of the expected range of 0.66–0.75 (Brown et al., 2004; Kooijman, 2010; Hatton et al., 2019). Second, the mortality rate parameter ζ₁ was selected to be slightly negative (−0.07), different from the initial distribution (mean 0.54). Note that the negative value for ζ₁ does not indicate a negative mortality, since ${\mathrm{\Lambda }}_k\propto e^{{\mathit{\zeta }}_{\mathrm{1}}}$ in Eq. (4). The mortality scaling (0.51) was also smaller than the mean value (0.54). To account for a large uncertainty in the trophic efficiency (α) (Eddy et al., 2020), we expanded its prior range to [0.06,0.4] (compared to previous estimates in Carozza et al., 2017). However, the optimization persistently selected values comparable to BOATSv1, with a mean of 0.14 (range 0.09–0.14; Table 3). Although separate pelagic and demersal communities could have different trophic efficiencies (Stock et al., 2017; Du Pontavice et al., 2020), for simplicity we adopt the same value here. The robustness of the optimized trophic efficiency suggests that sources of variability could be captured by other model parameterizations, e.g., the representative size of primary producers or the temperature dependence of growth and mortality. Lastly, growth activation energies (ω_a,A) for pelagic (0.50) and demersal (0.37) communities are larger than the prior values. Although the temperature dependence of mortality (ω_a,λ) is not significantly different from the initial values, the optimized values suggest a stronger sensitivity of growth compared to mortality for the pelagic community (${\mathit{\omega }}_{\mathrm{a},\mathrm{A}}-{\mathit{\omega }}_{\mathrm{a},\mathit{\lambda }}=+\mathrm{0.047}$ eV) and a stronger sensitivity of mortality for the demersal community (−0.082 eV).

Covariations between parameters in the 42-member optimized ensemble reveal compensations between parameter pairs (see Fig. 6a). The most significant compensations are between parameters controlling the biomass flow through the size spectrum and biomass losses (see Fig. 6b–d). For instance, an increase in the trophic efficiency (α) can be compensated for by a smaller predator–prey biomass ratio (β), which lengthens the food web (Fig. 6b; $r=-\mathrm{0.66}$). When more biomass flows across trophic levels, longer food chains ultimately lead to greater losses over the food web and thus similar fish biomass production. Alternatively, an increase in the trophic efficiency can be compensated for by an increase in the mortality parameter (ζ₁; Fig. 6c; r=0.42). Conversely, the mortality parameter ζ₁ decreases when the growth scaling exponent (b) increases (Fig. 6d; $r=-\mathrm{0.49}$), instead of decreasing, because of indirect impacts on the asymptotic size (m_∞).

Correlations between parameters that differ between pelagic and demersal food webs can also reveal trade-offs, particularly in how activation energies collectively affect the two communities (see Fig. 6e–g). For instance, an increased temperature sensitivity of growth for the pelagic community ${\mathit{\omega }}_{\mathrm{a},\mathrm{A}}^{\mathrm{P}}$ is matched by an increased sensitivity of growth for the demersal community ${\mathit{\omega }}_{\mathrm{a},\mathrm{A}}^{\mathrm{D}}$ (Fig. 6e; r=0.74) and a shift of the sensitivity of demersal mortality compared to demersal growth ${\mathit{\omega }}_{\mathrm{a},\mathrm{A}-\mathit{\lambda }}^{\mathrm{D}}$ (Fig. 6f; r=0.65). Another relationship between communities is observed for losses. As the temperature dependence of mortality for the pelagic community (${\mathit{\omega }}_{\mathrm{a},\mathit{\lambda }}^{\mathrm{P}}$) increases, increasing losses, there is a concurrent decrease in the representative size at the base of the benthic food chain (log ₁₀(m_β). This extends the food chain length, increasing losses in the demersal community (Fig. 6g; $r=-\mathrm{0.40}$).

While not exhaustive, this parameter analysis suggests trade-offs between biomass production and dissipation in pelagic and demersal communities (see Fig. 6a for further details).

Figure 6 Emergent relationships between model parameters. (a) Pairwise correlations between model parameters for the 42-member optimized ensemble. (b–g) Scatter plots illustrating the relationships between the most strongly correlated parameters: (b) trophic efficiency (α) and predator–prey mass ratio (β); (c) trophic efficiency (α) and mortality parameter (ζ₁); (d) growth scaling exponent (b) and mortality parameter (ζ₁); (e) growth activation energy for the pelagic community (${\mathit{\omega }}_{\mathrm{a},\mathrm{A}}^{\mathrm{P}}$) and demersal community (${\mathit{\omega }}_{\mathrm{a},\mathrm{A}}^{\mathrm{D}}$); (f) growth activation energy for the pelagic community (${\mathit{\omega }}_{\mathrm{a},\mathrm{A}}^{\mathrm{P}}$) and difference between growth and mortality activation energies (${\mathit{\omega }}_{\mathrm{a},\mathrm{A}-\mathit{\lambda }}^{\mathrm{D}}$) for the demersal community; (g) mortality activation energy for the pelagic community (${\mathit{\omega }}_{\mathrm{a},\mathit{\lambda }}^{\mathrm{P}}$) and representative size of organisms at the base of the demersal food web (m_β). In panel (a), circles indicate p values < 0.05 and stars p values < 0.01. In panels (b)–(g), the lines show linear regressions for the 42 parameter values of the ensemble; Pearson correlation r and p values are reported in each plot. In panels (b)–(g), red dots indicate the 5 final best parameter values for BOATSv2, while the black dots show the remaining 5 parameter values among the 10 best; gray dots indicate all other parameter values.

5.2 Features of the simulated catch

5.2.1 Global catch

Relative to BOATSv1, BOATSv2 corrects for the overestimate of high seas catch (see Fig. 5) while maintaining similar skill in reproducing historical variations of LME fish catch (see comparison of model ensemble means with SAU and WAT reconstructions; Fig. 4). BOATSv2 also shows improved skill in hindcasting the spatial evolution of catch (see Fig. 7) and the offshore and equatorward expansion of fisheries (Swartz et al., 2010; Guiet et al., 2020). In the 1950s (Fig. 7a–c), higher-latitude shelf regions and productive upwelling regions contributed the most to global catch. In the 1990s (Fig. 7d–f), while high latitudes still produced large catches, subtropical regions were also significantly exploited, especially in shallow regions. Productive high seas areas also supported significant fishing. The expanded representation of ecological processes accounts for most improvements in the high seas, while updating economic processes only yields minor improvements (see Appendix G).

Despite the closer fit to observations, model biases remain, in particular low offshore catches in the western equatorial Pacific and excessive catches in the northern and southern Atlantic (see Fig. 7d). These biases are not improved by the economic update and are likely related to ecological factors (see Fig. G1 panel b vs. c). However, it remains unclear if biases could also result from historical interactions between ecosystems and fishing effort or from changing environmental conditions. Processes not included in the model, such as habitat alteration by bottom-trawling gear, additional constraints on habitats such as dissolved oxygen (Deutsch et al., 2020), fish migrations and movement (Barrier et al., 2023; Guiet et al., 2022; Lehodey et al., 2008; Watson et al., 2015), or management and regulation, likely play a role in these biases. While the two observational catch reconstructions used to calibrate the model show differences, likely related to different approaches (e.g., using bathymetry or not; compare Fig. 7e, f), biases in simulated catches are apparent when comparing with either reconstruction.

Figure 8 shows the differences (residuals) between maximum simulated and observed peak catches in each LME. While there is an overall improvement from BOATSv1 to BOATSv2 (compare Fig. 8a, b to c, d), areas of overestimated (e.g., Indian Ocean) or underestimated catches (e.g., northwestern Pacific) are correlated between the two model versions, suggesting structural biases in the model. It is possible that accounting for features of coastal habitats such as coral reefs and mangrove forests could reduce regional biases, especially in Southeast Asia (Tittensor et al., 2010). Representation of biodiversity also remains crude, and additional functional types with life histories that differ from those of fish, such as cephalopods, could be considered (Denéchère et al., 2024). Finally, some larger predators that dive to feed on the deep scattering layer experience environmental conditions that differ from those at the surface (Nuno et al., 2022; Braun et al., 2023). Accounting for this effect could help reduce model biases.

Figure 7 Observed and simulated catch in the 1950s and 1990s, C_50s–90s (in log ₁₀; $\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). (a, d) BOATSv2 simulation forced with exponentially increasing catchability over time compared to (b, e) SAU and (c, f) WAT observations.

Figure 8 Residuals of observed vs. simulated maximum catch per LME, C_LMEmax $(C^{\mathrm{sim}}/C^{\mathrm{obs}})$. (a, b) BOATSv1 simulation compared to SAU and WAT observations. (c, d) BOATSv2 simulation compared to SAU and WAT observations.

5.2.2 Pelagic vs. demersal catch

Separate pelagic and demersal energy pathways allow simulation of higher taxonomic diversity. At the global peak of the 1990s, a large fraction of simulated demersal catch is derived from high latitudes (Fig. 9), in general agreement with observations (van Denderen et al., 2018). At lower latitudes, modeled demersal catches are as abundant as pelagic catches in shallow regions or near seamounts, also consistent with observations. Significant biases remain, however, such as in the North Atlantic, where the simulated demersal catch fraction is lower than observed, and the eastern tropical Pacific, where the demersal fraction is overestimated. The latter bias could reflect the parameterization of iron limitation, which reduces accumulation of pelagic biomass in the eastern topical Pacific, an HNLC region (see Appendix B).

Figure 9 Observed and simulated pelagic catch fraction $C_{\mathrm{90}\mathrm{s}}^{\mathrm{P}}/(C_{\mathrm{90}\mathrm{s}}^{\mathrm{P}}+C_{\mathrm{90}\mathrm{s}}^{\mathrm{D}})$ (in %). (a) BOATSv2-simulated fraction compared to (b) SAU observations. Very low catch levels are masked (i.e., below 10⁻³ $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$).

5.2.3 Deepening of the catch

The historical expansion of fisheries is associated with fishing in increasingly deep waters, i.e., a deepening of the catch (Morato et al., 2006; Watson and Morato, 2013). This can be attributed to the need to find new profitable fishing grounds beyond more accessible coastal regions, as well as improvements in fishing technology. The catch density per depth stratum from observational reconstructions reflects such expansion (see Fig. 10a and b for pelagic and demersal catch in SAU data). The deepening of demersal catches is consistent with increasingly deep fishing grounds, while the deepening of pelagic catches indicates an expansion of fishing effort towards deeper regions offshore.

In the model, decreasing biomass with depth slows the historical deepening of demersal catches (Fig. 10d). Similarly, higher costs and reduced catchability at greater depths delay the offshore expansion and deepening of pelagic catches (Fig. 10c). These factors collectively contribute to the slower development of fisheries in deep waters and the reduced catch fraction from the high seas, consistent with observational reconstructions (compare with BOATSv1 in Appendix H).

Figure 10 Fishing over increasingly deep seafloors. (a, b) Observed mean catch density (in log ₁₀; $\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) over depth strata for pelagic (a) and demersal (b) catch from SAU. (c, d) Simulated mean catch density over depth strata for pelagic (c) and demersal (d) catch in BOATSv2 (see Appendix H for comparison with BOATSv1).

5.3 Features of simulated fishing effort

The modeled nominal fishing effort aggregated across the 55 LMEs broadly matches observations (see observations in Fig. 4b), with a slightly earlier decline that falls within the uncertainty range. This could indicate that the model's effort responds to biomass depletion faster than observed or that the model underestimates the resilience of exploited stocks. It could also reflect the lack of management and subsidies in the model, which influence profitability and the progression of fishing effort.

The significant correlation between modeled and observed effort at peak catch across LMEs (see Fig. 11a; r_LME90s=0.57) lends support to the model's assumption of open-access dynamics. However, significant deviations remain (Fig. 11b). For instance, the model overestimates effort in highly productive shelf regions near the mouth of major rivers such as the Patagonian Shelf or North Brazil Shelf, suggesting development of fisheries that is too rapid compared to neighboring regions. Biomass redistribution by currents, or fish stock migrations, could correct this bias. At larger scales, the model underestimates effort across Southeast Asia, consistent with the underestimated peak catch. However, effort around Australia is also lower than observed, while the model overestimates catches there (see Fig. 8). This mismatch suggests regional differences in economic drivers or missing key habitats such as mangroves or reefs. Efficient management might also play a role, although it is unlikely to be the sole driver across the entire region.

Figure 11 Observed and simulated nominal effort in the 1990s, E_90s. (a) Simulated effort density per LME compared to observations, E_LME90s. (b) Ratio between modeled and observed effort by LME at the global catch peak in the 1990s. In panel (a), red dots and lines show BOATSv2 output, and blue dots and lines show BOATSv1 output; numbers indicate the LMEs (see Appendix E).

5.4 Features of simulated biomass

5.4.1 Global biomass

In the absence of fishing, BOATSv2 estimates a commercial fish biomass of 1.9 Gt aggregated over LMEs, slightly larger than previous estimates from BOATSv1 (1.6 Gt in Bianchi et al., 2021). However, because the biomass in the high seas is lower, the biomass in LMEs accounts for 68 % of the global biomass (2.8 Gt). This is significantly larger than the 50 % of BOATSv1 (3.3 Gt in Bianchi et al., 2021). Thus, BOATSv2 suggests a 10 %–15 % smaller pristine biomass than BOATSv1. When fishing is included and forced by the historical catchability increase, the BOATSv2 commercial fish biomass aggregated across LMEs declines by about 50 % from 1950 to 2000. This is consistent with both BOATSv1 simulations and global observational estimates (Fig. 4d; also compare with Worm and Branch, 2012). Interestingly, at peak catch, the LME–HSE difference between model versions is compensated for by differences in fishing effort so that both LMEs and HSEs hold approximately 50 % of the global biomass at this point (0.9 Gt within LMEs and 0.8 Gt in HSEs with BOATSv2, similar to 0.6 and 0.5 Gt, respectively, with BOATSv1; Bianchi et al., 2021). The optimized parameters of BOATSv2 suggest that the global fish biomass is about 40 % pelagic and 60 % demersal, a partitioning which could be relevant for the biogeochemical cycling and carbon export effects of fish. Comparing BOATSv1 and BOATSv2, the similar relative biomass distribution at peak harvest and the similar magnitude of pelagic biomass would suggest comparable estimates of export and sequestration by sinking fecal pellets (Bianchi et al., 2021). However, further analyses is needed to differentiate the roles of pelagic and demersal communities and their historical depletion in carbon and nutrient sequestration (Cavan and Hill, 2022).

5.4.2 Regional biomass distributions

In the model, shallow shelves and upwelling systems sustain on average 3 times more biomass per unit area than the high seas (10. vs. 2.9 g m⁻² within and outside the LMEs, respectively; see also Fig. 12a for local biomass gradients). Validating these predictions is challenging because of observational limitations; however, recent compilations provide a new way to assess them.

For demersal fish, scientific trawl compilations are now available from the Northern Hemisphere at locations ranging from Alaska to Europe (Maureaud et al., 2024; van Denderen et al., 2023). Figure 12 shows that BOATSv2 accurately simulates the average biomass across these LMEs and captures the biomass increase from the Gulf of Mexico–Florida (GM–FL) to Europe (EU) and the North American west coast (NA–W) (see circles Fig. 12b). However, the model underestimates the observed range of variability: while observations vary over almost 2 orders of magnitude, simulated biomasses vary only over 1 order of magnitude. The model also overestimates biomass for the North American east coast (NA–E).

These biases might reflect temporal offsets in the depletion of fish biomass over time due to exploitation (see also Fig. 4d), with the model failing to capture relative differences across LMEs. Indeed, regions where the model overestimates fish biomass still have relatively high simulated rates of biomass decline in the 2000s (e.g., LMEs 1, 5, 12). This can eventually deplete biomass to the observed levels. Conversely, regions where the model underestimates fish biomass (e.g., LMEs 10, 14, 60) are areas in which simulated fishing effort caused an early biomass decline (see Appendix I). These temporal mismatches could reflect regional differences in the rate of development of fisheries that are not captured by the simple globally homogeneous exponential increase in technology or by the open-access assumption and lack of management. For example, effective management in Alaskan fisheries has prevented the phase of overfishing that has been common in industrial fisheries worldwide (Worm et al., 2009). Alternatively, considering biodiversity could help explain these differences. For instance, in NA–W, the dominance of semi-pelagic Alaskan pollock may lead to an underestimation of our exclusively demersal biomass. Conversely, in NA–E, shifts from demersal to pelagic communities due to fishing can explain the overestimation of demersal biomass (Choi et al., 2004). Our approach does not capture these interactions between pelagic and demersal communities. Finally, discrepancies may be exacerbated because observations cover only a portion of each LME (compare mean biomass densities at LME level vs. grid cells where simulations overlap with observations in Appendix I).

Aggregated biomass observations for pelagic stocks are scarcer than for demersal stocks. We compare model output with scientific trawl data for coastal pelagic species in the California Current and Gulf of Alaska (shown by squares in Fig. 12b) (Zwolinski et al., 2012). The model simulates overall higher biomass densities than observed, showing a wider ranges of values. Similar to observations, simulated pelagic biomass densities are lower compared to demersal biomass. A caveat to this comparison is that estimates of pelagic biomass remain significantly uncertain due to challenges in sampling the three-dimensional oceanic environment, variability and aggregation in fish populations, uncertainty in sampled depth ranges, net avoidance by pelagic fish, and the limited selectivity of pelagic trawls (Kaartvedt et al., 2012; Zwolinski et al., 2012). In addition, fish migrations can redistribute fish biomass across life stages in ways not captured by the model.

Figure 12 Observed and simulated biomass density. (a) BOATSv2 global mean biomass density in the 1990s, B_90s. (b) Demersal biomass density from bottom trawl surveys (circles) versus simulations and pelagic biomass density from surface coastal pelagic species trawl data (squares) versus simulations. Biomass densities are averaged across LMEs for the reference decade of the 2000s. In panel (b), colors indicate the four different regions shown by the LME boundaries in panel (a) and numbers the specific LMEs (see Appendix E). In panel (b), the dots indicate median values, and horizontal or vertical lines the 25–75th percentile range; the size of the dots indicates the relative size of the surface area sampled. Appendix I provides a further comparison of biomass time series by LME.

5.5 Implications of the model update

The inclusion of distinct energy pathways and spatially variable economic drivers in BOATSv2 has a limited impact on the evolution of coastal fisheries over time but has a large impact on simulated high seas fisheries. All else being equal, BOATSv1 and BOATSv2 ensembles show very similar LME-level catch at the global peak (see $C_{\mathrm{LME}\mathrm{90}\mathrm{s}}^{\text{SAU-WAT}}$ Fig. 5c, d) and comparable progression across LMEs from 1950 through 2000. Both are in good agreement with observations (Fig. 13). The key improvement of BOATSv2 is the representation of high seas fisheries, where catches are delayed and greatly reduced during the historical period, bringing the model closer to observations (see dashed lines in Fig. 13). However, as fisheries keep developing, BOATSv2 still overestimates fishing in the high seas (compare the dashed red line with observations in Fig. 13 or the increasing trend in simulations in Fig. 5a). This discrepancy suggests either an improper representation of the historical rate of catchability increase in the simulations or missing mechanisms, such as horizontal migrations that redistribute biomass from the high seas to the coast.

Finally, the separation between pelagic and demersal communities influences fish production rates because these communities respond to different environmental drivers (van Denderen et al., 2021; Fredston et al., 2023). Compared to BOATSv1, this change could influence the resilience of fisheries to fishing and/or climate change. It could also alter the response to regulation, although we anticipate similar dynamics as in previous work (Scherrer and Galbraith, 2020; Scherrer et al., 2020). A separation of pelagic and demersal energy pathways is likely to impact the effects of fish on biogeochemistry (Bianchi et al., 2021; Le Mézo et al., 2022).

Figure 13 Time series of catch in the global ocean. Historical annual catch in the high seas (HS) vs. coastal seas (CS). Gray lines show observational reconstructions from Pauly et al. (2020) (SAU, dark gray) and Watson (2017) (WAT, light gray). Blue lines show output from BOATSv1 and red lines from BOATSv2, forced only with exponentially increasing catchability over time. Coastal (solid lines) and high seas (dotted lines) catches are shown separately.

6 Conclusions

We introduce BOATSv2, an expanded version of the BOATS model that includes multiple added features, including separation of demersal from pelagic communities, and improves simulation of high seas fisheries. New model features have a limited impact in coastal regions so that BOATSv2 simulates dynamics and variability in catch and biomass over LMEs that are similar to BOATSv1. The expanded representation of functional and taxonomic diversity allows more detailed comparisons with observations. In some cases, this reveals new model biases, such as in the simulation of demersal catches and biomass in the western North Atlantic.

We attribute improvements in the simulation of high seas fisheries to the separation of pelagic and demersal energy pathways, supporting the importance of distinguishing these communities (Blanchard et al., 2012; Petrik et al., 2019; Du Pontavice et al., 2020). We also introduced parameterizations of spatially heterogeneous economic drivers, i.e., fishing costs and catchability, which further improves the match to observations in the high seas. However, choosing between different formulations for these drivers (i.e., depending on distance from the coast, d_coast, or depth of the seafloor, z_bot) was only possible by testing plausible functional forms and retaining those leading to the largest improvements against empirical data. While this selection was not exhaustive, our final formulation is consistent with a variety of new observational constraints, such as the increase in fishing effort over seamounts (Kerry et al., 2022) and the historical deepening of fishing as technology progresses (Watson and Morato, 2013). We acknowledge that the specific choice of functional forms for these parameterizations is not well constrained and will likely require future refinement against observational diagnostics.

Because of the more accurate representation of high seas fisheries in BOATSv2 relative to BOATSv1, the fraction of catch that takes place in the High Seas at the time of the global catch peak is reduced from 31 % to 11 %, bringing it closer to the observed 8 %–9 %. Similarly, the mean depth of the catch declines from 1698 m in BOATSv1 to 420 m in BOATSv2, aligning it more closely with the empirical estimate of 154–372 m. This update should help reduce model uncertainties in future projections (Galbraith et al., 2017; Lotze et al., 2019; Tittensor et al., 2021) and provide a more accurate representation of the role of fish in global biogeochemical cycles (Bianchi et al., 2021; Le Mézo et al., 2022). Future model improvements could include a representation of the migration of fish stocks (Sumaila et al., 2015), the role of diverse coastal environments such as mangroves, reefs, and lagoons (Tittensor et al., 2010), and a representation of distinct mesopelagic communities (Irigoien et al., 2014; St. John et al., 2016; Hidalgo and Browman, 2019).

Appendix A: BOATSv2 governing equations and parameters

Table A1 List of ecological and economic parameters in BOATSv2 for pelagic (P), demersal (D), or both communities. Parameters selected by the Monte Carlo procedure are reported by the range of best selected values. See more details in Carozza et al. (2016, 2017).

^* Estimation of the fraction of large phytoplankton production following Dunne et al. (2005).

Table A2 Variables and governing equations for the ecological module of BOATSv2 for pelagic (P), demersal (D), or both communities. See more details in Carozza et al. (2016, 2017).

¹ Estimated from net primary production and surface temperature (in °C) following Dunne et al. (2005). ² Correction of trophic scaling when reduced growth in iron-limited regions is activated.

Table A3 Variables and governing equations for the economic module of BOATSv2 for pelagic (P), demersal (D), or both communities. See more details in Carozza et al. (2016, 2017).

¹ Not detailed in the present model description; see Scherrer and Galbraith (2020). ² Correction of catchability or cost when spatial economic parameterization is activated.

Appendix B: Global variability in surface nitrate

Figure B1 Minimum monthly sea surface nitrate concentration (in µM) from the World Ocean Atlas (Locarnini et al., 2006).

Appendix C: Spatial variation of costs

The cost of fishing varies by fishing gear and by targeted fish community (Lam et al., 2011). To best constrain spatially variable costs we use estimates of these separate fishing costs in the high seas for the main gear types (98 % of total effort) following data reported by Sala et al. (2018). Table C1 summarizes these estimated costs. These compare with the BOATS default fishing cost of 5.85 $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Carozza et al., 2017; Galbraith et al., 2017). We defined spatially variable costs as a function of distance to the coast d_coast and depth of the seafloor z_bot. Figure C1a illustrates the profile of distance-dependent costs and Fig. C1b the profile of depth-dependent costs (in $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$).

Table C1 Cost of fishing the high seas based on estimates from Sala et al. (2018) for the year 2016.

Figure C1 Cost per unit effort profiles $c_k^{\mathrm{coor}}$ in the global ocean (in $\mathrm{USD}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) (a) as a function of distance to the nearest coast $c_k^{\mathrm{coor}}=c_k^{\mathrm{dist}}$ and (b) as a function of depth of the seafloor $c_k^{\mathrm{coor}}=c_k^{\mathrm{depth}}$.

Appendix D: Spatial variation of catchability

Technology coefficients vary between gears (Palomares and Pauly, 2019), and gears are predominantly used in separate regions of the global ocean (Kroodsma et al., 2018; Kerry et al., 2022). Ultimately this can lead to spatially heterogeneous catchability of fish resources.

In order to better constrain the catchability, we use the reported difference of technology coefficients by gear estimated in 1995 and a coarse estimation of the contribution of each gear to the global fishing effort from 2015 through 2020 as reconstructed by Global Fishing Watch (GFW; see Kroodsma et al., 2018). Depending on the functional type predominantly targeted by a gear, pelagic vs. demersal, we estimate the mean technology coefficients for pelagic species to be Of(P)=1.3 compared to Of(D)=0.9 for demersal species (see Table D1).

Based on the observation that a dominant part of fishing effort on pelagic species by longliners occurs near seamounts (Kerry et al., 2022), we adjust the spatial catchability as a function of the depth of the seafloor such that it varies from a minimum of Pr(z_bot)=0.8 over deep seafloors (e.g., for tuna seiners; Table D1) to Pr(z_bot)=2.4 in shallow regions (such that the global mean is 1.3). Figure D1 illustrates the reference profile of technology coefficients Pr(z_bot) used for the analysis.

Table D1 Technology coefficient per fish community. The coefficients per gear are based on reported values in Palomares and Pauly (2019). Each gear is linked to the dominant resource it targets, pelagic (Pel) or demersal (Dem), and associated with the fraction of global fishing effort from 2015 through 2020, as reported by Global Fishing Watch. We reported the mean technological coefficient weighted with effort by gear, when available.

Figure D1 Relative technology coefficient profiles Pr(z_bot) in the global ocean.

Appendix E: Large marine ecosystems and high seas ecosystems

Figure E1 Regional masks used to compare observations and simulations. (a) Large marine ecosystems. (b) High seas ecosystems adapted from Weber et al. (2016).

Appendix F: Pelagic and demersal catches

Table F1 The Sea Around Us (Pauly et al., 2020) functional types associated with pelagic and demersal catches.

Appendix G: Global catch distribution between BOATS versions

Figure G1 Simulated catch in the 1990s, C_90s (in log ₁₀; $\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). (a) BOATSv1. (b) Updated version with improved ecology – BOATSv2-Bio. (c) Final update including improved economics – BOATSv2.

Appendix H: Historical catch deepening

Figure H1 Fishing over increasingly deep seafloors. (a) Observed and (b, c) simulated mean total catch density (in log ₁₀; $\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) over depth strata. Compared to observations (a), BOATSv1 (b) fails to capture the deepening, while BOATSv2 corrects it.

Appendix I: Historical biomass variation in selected LMEs

Figure I1 Observed and simulated biomass decline in LMEs for demersal and pelagic biomass. Each panel shows the simulated historical fish biomass density (in g m⁻²) decline, averaged across the selected LME (dotted line) or averaged over 1° grid cells where observations are available (plain line). These are compared to the range of observed biomass density per LME over years in the 2000s, indicated by the median value (plain black lines) and the 25th and 75th percentiles (dotted black lines). Colors indicate neighboring LMEs, North American LMEs along the east coast (green), North American LMEs along the west coast (light blue), the Gulf of Mexico and Florida LMEs (in orange), and European LMEs (dark blue).