Modeling the Population Effects of Hypoxia on Atlantic Croaker (Micropogonias undulatus) in the Northwestern Gulf of Mexico: Part 2—Realistic Hypoxia and Eutrophication

Quantifying the population-level effects of hypoxia on coastal fish species has been challenging. In the companion paper (part 1), we described an individual-based population model (IBM) for Atlantic croaker in the northwestern Gulf of Mexico (NWGOM) designed to quantify the long-term population responses to low dissolved oxygen (DO) concentrations during the summer. Here in part 2, we replace the idealized hypoxia conditions with realistic DO concentrations generated from a 3-dimensional water quality model. Three years were used and randomly arranged into a time series based on the historical occurrence of mild, intermediate, and severe hypoxia year types. We also used another water quality model to generate multipliers of the chlorophyll concentrations to reflect that croaker food can be correlated to the severity of hypoxia. Simulations used 100 years under normoxia and hypoxia conditions to examine croaker population responses to the following: (1) hypoxia with food uncoupled and coupled to the severity of hypoxia, (2) hypoxia reducing benthos due to direct mortality, (3) how much hypoxia would need to be reduced to offset decreased croaker food expected under 25 and 50% reduction in nutrient loadings, and (4) key assumptions about avoidance movement. Direct mortality on benthos had no effect on long-term simulated croaker abundance, and the effect of hypoxia (about a 25% reduction in abundance) was consistent whether chlorophyll (food) varied with hypoxia or not. Reductions in hypoxia needed with a 25% reduction in nutrient loadings to result in minimal loss of croaker is feasible, and the croaker population will likely do as well as possible (approach abundance under normoxia) under the 50% reduction in nutrient loadings. We conclude with a discussion of why we consider our simulation-based estimates of hypoxia causing a 25% reduction the long-term population abundance of croaker in the NWGOM to be realistic and robust.


Introduction
Quantifying the population-level effects of hypoxia on coastal fish species has been challenging. Fish populations exhibit complex dynamics that are not easily attributed to a single factor (Rose 2000). Many coastal fish populations live tens of years and their life stages use different habitats. This results in individuals being affected by a multitude of environmental and Communicated by Dennis Swaney biological factors, whose influence can vary within years (e.g., seasonally) and among years (interannual). Quantifying population-level effects of hypoxia is further complicated by the stressor itself (low dissolved oxygen (DO)) varying spatially and temporally, and combined with fish movement patterns, resulting in difficulties in determining time-dependent exposures of individuals in the field. Also, while very low DO undoubtedly results in mortality of individuals, most fish attempt to avoid lethal DO levels, and thus, the realized effects of low DO depend on realistic modeling of behavior and dealing with sublethal (e.g., growth) and indirect effects (e.g., displacement to poorer habitat). This same suite of issues creates challenges in quantifying population responses to hypoxia whether the analysis is based on long-term monitoring data or modeling ).
Our population modeling strategy is to focus on a wellstudied species that was known, based on its observed avoidance behavior, to experience exposure to low DO conditions that caused direct growth and reproduction effects in the laboratory. We then built up the population model, using a spatially explicit individual-based approach, to allow for simulation of realistic exposures and to use laboratory information to clearly link the time-varying exposures to growth, mortality, and reproduction effects on individuals. We described the individual-based population model in Rose et al. (2017a) and the exposure-effects model specifically designed to deal with fluctuating DO exposures in Neilan and Rose (2014). As part of the Rose et al. (2017a) analysis, we used idealized (climatological) environmental conditions repeated yearafter-year so that long-term population responses to hypoxia could be easily identified and the reasons for the responses within the model understood. In this companion paper, we focus on expanding that analysis to include more realistic spatial and temporal variation in DO concentrations and to include possible changes in croaker food (river-related production; mortality of benthos) that can occur concomitantly with the severity of hypoxia.
Low DO concentrations associated with the summer hypoxia in the northwestern Gulf of Mexico (NWGOM) exhibits considerable spatial and temporal variability. The midsummer areal extent of hypoxia during 1985 to 2011 varied from 700 to 23,200 km 2 , Table S2). The location of the hypoxic zone, based on the same snapshot spatial cruises in late July, typically not only included a core area in many years but also showed wide variation in the location of the edges of the core zone and in the number and locations of patches of higher DO (10's kilometers scale) within the general hypoxia zone (Rabalais et al. 2007DiMarco et al. 2010). Based on long-term monitoring at stations located in the center of the hypoxic zone (C transect), hypoxia has occurred starting as early as late February and lasting as late as October, with most years showing persistent hypoxia between June and September (Rabalais et al. 2007;Babin 2012). Low DO also shows finer-scale temporal and spatial variability that is superimposed upon the seasonal-scale variation. DO concentrations on the shelf monitored at key locations have been observed to change on the scale of hours (Babin and Rabalais 2009;Bianchi et al. 2010a), and cross-shelf variability in the distribution of hypoxia creates mesoscale (km level) patches of hypoxia (DiMarco et al. 2010). The 15 min to 1 h monitoring of near bottom DO concentrations within the hypoxic zone showed rapid increases in DO due to various types of mixing events (Babin 2012) and other local hydrodynamic influences (Bianchi et al. 2010a). Using individual samples, vertical thickness of the hypoxic zone varied from less than a meter to 20 m over the historical record ( Fig. S2 in Obenour et al. 2013). Near continuous measurements of DO using a towed body (Scanfish) that undulated between 2 m below the surface and 2 m above the bottom showed that changes in DO (e.g., 1-3 mg l −1 ) occurred on the scale of a kilometers (Roman et al. 2012;Zhang et al. 2014). We did an analysis of the same Scanfish data to further qualify meter-scale variability in which we only looked at a specific 1-h time series (measurements every 0.5 s) when all values were less than 4 mg l −1 and within 5 m of the bottom. We sampled these values (1100 points) at different distances from each other: 25% of the values showed a change in DO of 0.45 mg l −1 or higher over 20-m distance. These multiple scales of temporal and spatial variability in DO underlie the relatively smooth spatial maps of hypoxia (e.g., Obenour et al. 2013) that serve as annual indicators of the location and degree of hypoxia.
Understanding the link between degree of hypoxia and possible food availability to croaker is complicated. We focus on two possible linkages: (1) low DO causing direct mortality of benthos and (2) food being linked to hypoxia conditions via river flow-driven productivity whereby high nutrient fluxes via riverine loadings lead to more food and more severe hypoxia. Hypoxia effects on benthos (e.g., mortality, body size) has been documented for the NWGOM Qu et al. 2015) and other coastal systems (e.g., Levin et al. 2009;Kodama and Horiguchi 2011).
The evidence for the correlation between hypoxia, primary production, and river flow is quite strong. Hypoxia extent has been related to nitrogen loadings as a threshold response in several well-studied systems (Conley et al. 2009) and Hondorp et al. (2010) found weak but consistent relationships between nutrient loadings and severity of hypoxia across 23 ecosystems. Nitrogen loadings from the Mississippi River have been correlated to the area of the hypoxic zone in the NWGOM (Turner et al. 2008), and used in regression models to predict hypoxic area (Turner et al. 2012). Nutrient loadings from the river is also correlated with the primary production associated with the plume in the NWGOM (Lohrenz et al. 1999;Rabalais et al. 2002;Dagg and Breed 2003;Wysocki et al. 2006), as well across the broader Louisiana and Texas shelf (Rabalais et al. 2007) and in other river-influenced marine systems (Rabalais 2002). Remote sensing data for 2000 to 2003 showed interannual and monthly relationships among river flow, nutrient loadings, and chlorophyll concentrations over large area of the NWGOM; however, there was not a simple relationship detected between river discharge or nitrate loading with DO concentrations measured at a station within the hypoxic zone (Walker and Rabalais 2006). The next logical step from primary production to zooplankton is generally supported by the empirical information. Zooplankton (both micro and copepods) consume a large portion of the phytoplankton and so it would seem logical that higher primary productivity would lead to higher zooplankton; indeed, copepod densities are relatively high in the plume-influenced waters (Dagg and Breed 2003). Zhang et al. (2014) showed a rough correspondence between areas of high chlorophyll concentrations and high zooplankton densities on summer time transects in the NWGOM. Subsequent analyses showed a relatively strong rank correlation (6 of 12 correlations were greater than 0.5) between measured chlorophyll concentrations and zooplankton densities with values paired along these multiple transects (H. Zhang, personal communication). Peterson et al. (2002) documented a similar correspondence of high zooplankton fecundity rates and chlorophyll concentrations in coastal Oregon.
The final step of primary production and zooplankton subsequently affecting the benthos is more tenuous. Rabalais (2002) proposed a conceptual model that showed benthos productivity (as well as demersal fish yield) increasing with nutrient loadings until finally dropping at high nutrient loadings. Herman et al. (1999) found a positive relationship between primary production and macrobenthos biomass across different estuarine systems. Cloern (2001) used conceptual models and basic ecological food chain theory based on a review of the features and productivity of many systems and concluded that there was clear evidence for nutrient enrichment affecting phytoplankton biomass and primary production and also evidence in some systems for effects to further influence the benthos. Marcus and Boero (1998) not only cited several studies documenting pelagic fueling of benthos production but also added how the resuspension of nutrients from the benthos (rather than river-derived nutrients) can also act as the trigger for the higher pelagic production that then sinks to fuel the benthos. Empirically based qualitative and correlative relationships between pelagic production and benthos production have been shown for generally colder and deeper systems (e.g., Ambrose and Renaud 1995;Grebmeier et al. 1988;Johnson et al. 2007;Gooday 2002). Rowe (1971) and Rowe et al. (1974) analyzed data across several large coastal systems and determined that depth was most important in determining deep water benthos biomass but that surface productivity was also an important factor. In a review of global patterns, Wei et al. (2010) concluded that surface primary production was an important driver of benthic biomass on the sea floor.
Another piece of evidence, albeit circumstantial, for linking hypoxia and croaker food via nutrient loadings is the crosssystem comparisons of nutrient loadings versus fishery landings (e.g., Nixon and Buckley 2002;Caddy 2000). These comparisons show a clear relationship between nutrient loadings to an estuary and various indices (often catch or landings) of fish species. Breitburg et al. (2009) looked across about 24 different coastal systems and showed either increasing landings or a bowshaped relationship over the observed range of annual nitrogen loadings. This included fish species divided into pelagic planktivores, bentho-pelagics, and benthics. The key for our analysis here is whether these relationships also apply to how variation in annual nitrogen loadings would affect interannual differences in fish biomass within a given system.
In this paper, we use the individual-based croaker population model described in Rose et al. (2017a) and simulate long-term population responses under realistic variation in DO and with croaker food coupled to the severity of the hypoxia. We use DO concentrations from a 3-dimensional (3-D) coupled hydrodynamics-water quality model (Justić and Wang 2014) for the NWGOM as input to the croaker model. We also vary the food availability to croaker in two ways. The first is hypoxia reducing benthos due to direct mortality. The second way assumes food is positively correlated to the severity of hypoxia because both increase with increasing riverine nutrient loadings. We use the results of a second 3-D water quality model (Fennel et al. 2011(Fennel et al. , 2013 to generate how chlorophyll (assumed an indicator of croaker food) would respond to reduced riverinebased nutrient loadings. These changes were then imposed on food availability of croaker in the population model based on the severity of hypoxia each year in simulations. Reduced nutrient loadings therefore have competing effects on croaker: less nutrients reduce food and therefore reduce croaker while the same reduced nutrients act to lessen hypoxia that increases croaker. We simulated reduced nutrient loadings (less croaker) and then determined how much hypoxia would need to decrease (more croaker) to offset the food-related decrease in croaker. Finally, we performed a sensitivity analysis in which we changed key assumptions about avoidance movement, and isolated growth, mortality, and reproduction effects of hypoxia. We conclude with a discussion of why we consider our simulation-based estimates of the long-term effects of hypoxia on croaker both realistic and robust, and the next steps in our continuing quest to add more realism to model simulations.

Model Description Overview and Grid
The individual-based population model simulates the hourly growth, mortality, reproduction, and movement of individual croaker on a 2-D spatial grid. The grid consisted of 1 km × 1 km cells within a 800 × 300 rectangle that covered the area between Mobile Bay, AL and East Matagorda Bay, TX, and extended from about 70 to 200 km offshore (Fig. 1). The grid was centered on Louisiana because croaker occurs in high biomass there (Monk et al. 2015) and to capture the dominant summertime hypoxia event often located on the Louisiana shelf. A model year was 365 days long and began on September 1. The model is described in detail in Rose et al. (2017a).
Hourly temperature and surface chlorophyll-a concentrations were specified for each cell for a year and that year's set of hourly values were repeated every year (see ESM-1 in Rose et al. 2017a). Temperature was specified from estimated climatological values and we used bottom values for cells until 100-m depth was reached, after which the values for 100-m depth were used for cells. Temperature affected growth and movement, and chlorophyll-a concentration was used as an index of food availability. Temperature and chlorophyll-a started with mean monthly values, and these were interpolated to the IBM grid. Temperature was interpolated to daily values, which were used for all hours within each day, while monthly chlorophyll-a values were used for all hours and days within a month. Each cell was also assigned a DO concentration, to simulate different types (location, duration, patchiness) of hypoxic conditions. In this paper, we used the output of DO concentrations from a 3-D hydrodynamics-water quality model (Justić and Wang 2014) to specify the hourly DO concentrations on the IBM grid during the summer months.

Super-Individuals
The model used a super-individual approach (Scheffer et al. 1995;ESM-3 of Rose et al. 2017a). Each super-individual is assigned an initial worth when it is added to the population, which is the number of identical population individuals it represents. Mortality reduced the worth of over time, and super-individuals were removed from the model at age 8. All computations and output variables related to populationlevel metrics such as abundance, biomass, weight, fecundity, and spatial distributions are based on an individual's attributes multiplied by its worth. The model is described using the growth, mortality, reproduction, and movement of an individual; all computations in the model account for the worth of the super-individuals.

Individual Processes
Each individual croaker was tracked hourly as it progressed through the life stages of egg, yolk-sac larva, ocean larva, estuary larva, early juvenile, late juvenile, and adult. Superindividuals were started on each day there was spawning. Each individual was tracked (initially as an egg) as it experienced grid-wide temperatures that determined on what hour it switched from egg to yolk-sac larvae, yolk-sac larva to ocean larva, and ocean larva to estuary larva. When an individual was ready to go from estuary larva to early juvenile, it was placed on the 2-D grid with a starting length of 32 mm and a worth based on the mortality rate experienced from starting as a new super-individual representing eggs until the end of the estuary larval stage. Development rates were temperature-dependent, while mortality rates were stage-dependent. Early juveniles were restricted to cells shallower than 5 m.
Individuals on the grid were followed hourly through growth, mortality, reproduction, and movement. Growth, mortality, reproduction, and movement depended on the local conditions in their cell, and in some cases (e.g., avoidance movement), on conditions in surrounding cells.
Growth each hour was a bioenergetics equation with chlorophyll-a in the cell determining realized consumption: where P-value is the fraction of maximum consumption (C max ) realized, R is respiration, SDA is specific dynamic action, Eg is egestion, Ex is excretion, and ρ p and ρ f are the energy densities (J g −1 ) of consumed prey and of croaker. Maximum consumption and basal respiration rate were computed based on allometric functions of weight and of temperature. Egestion (Eg) was computed as a fraction of consumption, and SDA and excretion (Ex) were computed as fractions of consumption minus egestion.  The P-value was computed each hour from the average daily chlorophyll-a concentration Chl in each cell (denoted by column i and row j): The dependence of P-value on chlorophyll concentration ( Fig. 2) was designed so that croaker maintain reasonable consumption rates over a broad range of concentrations, likely due to their behavioral adjustments to prey availability, but then consumption rises rapidly (patch of dense prey) or drops quickly (diluted prey densities) outside of this broad range.
The daily growth from Eq. 2 was then divided by 24 and used to update weight (and then length) of the individual for that hour. Early juveniles became late juveniles when they exceeded 97.5 mm, and late juveniles became adults at 180 mm. Values for the bioenergetics parameters were specified for three size ranges of croaker (see Table 2 and ESM-4 in Rose et al. 2017a).
Mortality rate was specified as a stage-specific constant background rate, starvation, old age, and from exposure to low DO. Different background rates for estuary larva through adult stages were specific to east and west of 94°longitude. Late juvenile mortality rate was specified as density-dependent, with mortality rate increasing with increasing abundance of late juveniles in the same cell.
On September 1 of each year, the maturity of each individual was determined based on its length, and if mature, it was assigned a degree-days value that determined when they would start spawning in the fall and winter of the upcoming year. September 1 is after summer growth and before earliest spawning. When an individual accumulated its assigned degree days, it initiated spawning. Each individual released 12 batches of eggs: one batch every 3 days over a consecutive 34day period. Eggs per batch depended on weight of the individual at the time of spawning. Spawning was concentrated in simulations from October through March.
Individuals moved hourly based on a kinesis algorithm, and when they entered a cell with low DO, they switched to an avoidance algorithm. Kinesis used the sum of inertial and random velocities to determine the velocity of the individuals in the x and y directions on the grid. The position was then updated by multiplying the velocities by the 1 h time step and obtaining new x and y locations (and potentially a new cell identifier). Kinesis used temperature as the cue; temperature in a cell close to the specified optimal temperature for that stage resulted in the inertial velocity component being weighted more importantly in the sum of inertial and random velocities, while temperatures away from the optimal resulted in heavier weighting on the random velocity component.
Avoidance was triggered when an individual arrived at a cell determined at the beginning of the hour from kinesis had a DO concentration below a threshold (e.g., 2.0 mg l −1 ). Avoidance was initiated before the individual experienced the low DO value. The avoidance algorithm involved the individual assessment of the DO and temperatures in a neighborhood of cells centered on the cell determined from kinesis (i.e., that had low DO that triggered avoidance) and then moving towards the cell with DO greater than the threshold and with a temperature closest to their optimal value. Individuals were exposed to DO concentrations above the threshold (i.e., no avoidance triggered) that still cause growth and reproduction effects, and less frequently to DO values less than the threshold (cause mortality along with growth and reproduction effects) because avoidance was not completely effective every time. Under both kinesis and avoidance, early juveniles could only move to cells with depths less than 5 m, late juveniles to cells with depths less than 20 m, and adults to cells with depths greater than 5 m.
As individuals move with kinesis and avoidance, they experience the conditions in the cells and get exposed hourly to DO concentrations. We used three exposureeffects submodels that related hourly exposures to DO to reduced growth, increased mortality, and reduced annual fecundity (Neilan and Rose 2014). Each hour, the immediate effect of the DO for that hour was computed for growth (f G ), mortality (f M ), and reproduction (f R ) using the same general formulation: where DO(t) is the DO concentration in the cell, DO NE is the DO threshold below which an effect occurs, and α and β are non-negative constants. The values from Eq. 3 were used to determine the reductions in hourly growth, hourly survival, and, averaged over 10 weeks, to determine the reduction in all egg batches from that individual for the next spawning season. The DO exposure effects were only evaluated for late juveniles, age 1 adults, and age 2 adults, and all effects were reset each year so any accumulation of effects did not carry over to the next year. While Eq. 3 was evaluated every hour for the entire year, effects only occurred when DO concentrations were below the no-effects concentrations (i.e., summer months).

Simulated Dissolved Oxygen
A coupled hydrodynamics-water quality model (Wang and Justić 2009;Justić and Wang 2014) was used to generate DO concentrations as input for the IBM. The hydrodynamics model was an implementation of FVCOM; the WASP model was used to simulate water quality. The domain for FVCOM model extended from Mobile Bay, AL, to an area in the Gulf of Mexico south of the mouth of Galveston Bay, TX, and covered the area off the Louisiana-Texas Shelf affected by hypoxia (gray area in Fig. 1). Validation of the DO predictions from the FVCOM-WASP model was carried out using 2002 as the reference year (see Justić and Wang 2014), and then two additional years (2001 and 2010) were simulated to obtain 3 years of DO concentrations for use as input to the IBM.
FVCOM is a 3-D, primitive equation coastal ocean circulation model that features an unstructured triangular grid to resolve complex coastlines and a terrain-following coordinate system to convert irregular bottom topography into a regular computational domain (Chen et al. 2003(Chen et al. , 2007. The model computational grid consisted of 14,740 triangular nodes and 28,320 triangular elements, and had a variable 1-10 km horizontal resolution depending on shoreline proximity. The vertical grid domain consisted of 31 uniformly distributed σ layers and had a 0.2-2.0-m resolution within the general hypoxia area and vertical resolution of about 10 m in deeper waters on the shelf. The WASP water quality model simulated nitrogen and phosphorous kinetics, carbonaceous biological oxygen demand in the water column and sediment, phytoplankton biomass expressed as chlorophyll, and DO dynamics. Water quality variables were outputted hourly for each node of the hydrodynamics model grid.
We used the simulated hourly DO concentrations for each node for June 1 through August 31 of each of the 3 years. Each FVCOM-WASP simulation started on April 1, and the first 60 days were ignored for spin-up and because of little hypoxia in May. The rectangular IBM grid was then overlaid on the bottom layer of the FVCOM grid, and each IBM grid cell was assigned the DO value from the triangle that contained the center location of the IBM cell.
Simulation years 2001 and 2002 from FVCOM were considered years of severe hypoxia, and 2010 was considered an intermediate hypoxia year. All calculations about DO from the FVCOM presented here use the DO concentrations after they were interpolated to the IBM grid. Based on available monitoring data, Rose et al. (2017a) used definitions of hypoxia severity based on July area: mild (<10,000 km 2 ), intermediate (10,000-20,000 km 2 ), and severe (>20,000 km 2 ). With the FVCOM simulations of DO, we can examine daily values over the entire grid throughout the summer period, which provides additional information to determine the annual severity of hypoxia in addition to information just from a cruise during the the last 2 weeks in July. We therefore computed the average hypoxic area in two ways: maximum 10-day moving average during the summer months and the last 2 weeks in July. The values for the maximum 10-day average (and for late July) August (10,000-15,000 km 2 ) that would be considered intermediate that were not reflected in the July only average area (Fig. 3a). Even though the July hypoxia areas in years 2001 and 2002 were less than our previous definition of severe (i.e., >20,000 km 2 ), we considered then to be severe years because both showed a general increase in hypoxia area throughout the summer, with area hovering around 15,000 km 2 after July 20 and both had a week or more near 20,000 km 2 during August. Snapshots in mid-month of June, July, and August show the generally lower DO occurring over broader areas for our designated severe years compared to the intermediate year ( Fig. 4). We still used the earlier definitions of severity based only on late July areas to determine the frequency of mild, intermediate, and severe years because many years are needed to estimate the frequencies. The snapshots also demonstrated the kilometer-scale variability in the simulated DO concentrations. Time series of the DO concentrations used in the IBM at a key location (C6, triangle in Fig

Eutrophication and Croaker Food
We used multipliers applied to daily values of chlorophyll to simulate how riverine nutrient loadings could lead to a positive relationship between food availability to croaker and severity of hypoxia (i.e., more nutrients leads to more food and more severe hypoxia). The multipliers were derived from simulation results of a ROMS water quality model for the Gulf of Mexico (Fennel et al. 2011(Fennel et al. , 2013Laurent and Fennel 2014). They performed simulations for 2001-2007 under current conditions, and repeated the same years but under a 25% reduction and a 50% reduction in inorganic forms of nitrogen and phosphorus in the Mississippi and Atchafalaya Rivers (Laurent and Fennel 2014). Each of the ROMS years from 2001 to 2007 was categorized into mild, intermediate, or severe hypoxia based on their late July hypoxic areas and the same criteria as used to categorize FVCOM simulation years. We then grouped ROMS' years together based on their hypoxia type (mild: 2003, intermediate: 2004 and 2006, severe: 2001, 2002, 2005, and 2007). The year 2005 was difficult to easily categorize because while its late July area (11,800 km 2 ) was intermediate, more detailed examination of the DO dynamics showed periods of consistently low DO in early July and August. We therefore analyzed the ROMS chlorophyll output with 2005 treated as intermediate and alternatively as severe.
Because we used multipliers to adjust chlorophyll to allow for differences between the ROMS simulation of absolute chlorophyll concentrations and our climatological chlorophyll concentrations, we needed to decide which hypoxia year type was considered baseline for chlorophyll. We assumed our climatological chlorophyll in the IBM simulations was representative of mild hypoxia years. We examined specific years that comprised the climatological values, and qualitatively judged that our climatological chlorophyll was similar to mostly mild and some intermediate hypoxia years, and thus, we used the chlorophyll in the mild hypoxia ROMS year (2003) as the denominator for deriving the multipliers. The multipliers were not sensitive to whether we used mild or intermediate years as the denominator because the chlorophyll concentrations in the mild and intermediate years, once averaged for our use, were roughly similar.
We examined several schemes for deriving our multipliers from ROMS-generated average monthly values of surface layer chlorophyll. These included dividing the IBM grid into two subregions (inshore and offshore, Fig. 1) as well as six subregions (each offshore and inshore further subdivided into three boxes). For each scheme, we first averaged the monthly values for each year for spring (March through May) and for summer (June through August), and for each subregion.   After examination of the ratios by seasons using various spatial resolutions, we decided to use a simple set of multipliers. We used a single multiplier for all of the three inshore boxes (squares in Fig. 5). All offshore subregions (less influenced by the nutrient loadings) were reasonably set to a multiplier of 1.0. Under baseline conditions, we kept chlorophyll concentrations the same in mild years and increased chlorophyll-a in each cell each time step by 10% for intermediate and by 20% for severe hypoxia years. We reduced the multipliers the same percent for mild, intermediate, and severe years within each nutrient reduction and season combination: 5% for spring and 10% for summer for 25% nutrient reduction and 15 and 30% for the 50% nutrient reduction. The exception was the 17% reduction (1.2 to 1.0) rather than 10% for severe years in summer under the 25% reduction.
These derived multipliers used in simulations generally were within the range predicted by the output of the ROMS model (years 2001-2007 divided by 2003) except for the 1.2 multiplier for severe years in summer under baseline conditions (circles versus adjacent squares of the same color in Fig. 5). We used the larger increase in summer than suggested by the ROM results because other spatial averaging schemes and dividing by intermediate years resulted in multiple values that sometimes exceeded 1.1. Also, the high productivity of the Mississippi and Atchafalaya Rivers plumes would be highest in our inshore region (Lohrenz et al. 1997;Rabalais et al. 2002). The percent reductions are larger than the multiplier effect themselves because the multipliers for 25 and 50% predictions are compared to baseline that itself gets adjusted for hypoxia year type. For example, under the 50% reduction in nutrients and a severe hypoxia year, chlorophyll-a would be reduced by multiplying baseline values by 0.8; under baseline nutrients, a severe year would be multiplied by 1.2; therefore, the reduction is from 1.2 to 0.8 (or 33%).

Sequence of Hypoxia Year Types
Model simulations were for 140-years that used a random sequence of mild, intermediate, and severe hypoxia years for years 41 to 140. We used the same process to derive the probabilities of mild, intermediate, and severe hypoxia years as used in Rose et al. (2017a): 0.18 for mild, 0.52 for intermediate, and 0.30 for severe. In the IBM simulations, we assumed normoxia for mild years, 2010 from FVCOM for intermediate years, and years 2001 and 2002 from FVCOM were equally likely for severe years.
We used a single replicate of the sequencing of mild, intermediate, and severe hypoxia year types for all analyses in this paper. As in Rose et al. (2017a), we initially used 12 replicate 140-year time series that had different sequences of hypoxia year types to account for any effects on the results due to the different ordering of mild, intermediate, and severe hypoxia years (see Fig. 6 in Rose et al. 2017a). Using FVCOM-derived DO concentrations, there were little differences in averaged values of croaker abundance and other model outputs across the 12 replicates. For example, averaged croaker abundance (age 2 and older), computed over years 100 to 140 varied by less than 5% (7.051 versus 7.732 × 10 6 ) across the 12 replicates, and the CVof the 12 40-year annual abundances (within replicate variability) only varied from 2.7 to 5.7%. The variability among replicate simulations with coupled chlorophyll was similarly small. Average annual abundances of croaker for three replicates with hypoxia and chlorophyll coupled were 11.04, 10.99, and 10.67 million for no hypoxia effects and 8.08, 8.22, and 7.78 million with hypoxia effects. We used replicate 1 from Rose et al. (2017a), which had 22% mild years, 55% intermediate, and 23% severe years for all analyses here. For the 40 years (101 to 140) used in model output presentations, replicate 1 had 20% mild, 53% intermediate, and 27% severe years. Using one replicate helps simplify the presentation of model results without loss of generality.

Simulation Set 1: Isolating the Hypoxia Effect
The first set of simulations was designed to isolate the effects of hypoxia on the croaker population. Two conditions were simulated. The first condition was with chlorophyll uncoupled to hypoxia, and thus, the climatological chlorophyll concentrations (still varied within each year) were used every year. This allowed evaluation of hypoxia effects using FVCOM DO values when all other conditions are climatological. The second condition included chlorophyll being adjusted by the spring and summer inshore multipliers (Fig. 6) based on the severity of hypoxia. Chlorophyll was adjusted every year based on the hypoxia year type. Coupling food and hypoxia enables assessment of hypoxia effects when their effects can offset each other. For each of these two conditions, a simulation was performed without hypoxia (normoxia every year) and with hypoxia varying from year-to-year (replicate 1 pattern) using FVCOM output. Comparison of the without and with hypoxia for uncoupled and coupled conditions enables isolation of the hypoxia effect under both assumptions of chlorophyll being related or uncoupled to hypoxia.

Simulation Set 2: Hypoxia-Induced Benthic Mortality
The second set of simulations examined the population responses of croaker as result of hypoxia causing direct mortality on benthos (croaker food). Direct mortality of food was implemented by adjusting the P-value in each cell on each hour DO was less than 2 mg l −1 .
where P-value B is the P-value under normoxia, γ is the food recovery rate (P-value h −1 ), and Q is the hours required for 50% mortality of the food due to low DO. The P-value was then also constrained to be less than the baseline P-value (DO should not cause more food). A DO-adjusted P-value was computed for each cell over time based on the DO-adjusted P-value from the last hour and P-value from baseline (normoxia) for the next hour. On each hour, the smaller of the two P-values (DO-adjusted and baseline) were set to the P-value for that cell. This allows the effects of low DO to accumulate while also permitting the food effect to recover based on their own productivity. Q was set to 1000 (about 40 days), which was a representative of the times to 50% mortality for prey types eaten by croaker (Fig. 3 in Vaquer-Sunyer and Duarte 2008), and γ was set to 3.47 × 10 −4 based on production to biomass ratio value of 2 year −1 (Arrequin- Sanchez et al. 2002;de Mutsert et al. 2016) and assuming the production occurs for 240 days of the year. We also included an unrealistic simulation of 50% mortality occurring after only 4 days (Q = 100 h) to demonstrate that the model could predict a croaker population response. Based on mortality rates of benthos due to low DO, the 40-day assumption is the realistic scenario.

Simulation Set 3: Reduced Nutrient Loadings
The third set of simulations returned to the coupled chlorophyll-hypoxia conditions as in the first set of simulations, but now using the chlorophyll multipliers associated with the 25 or 50% reduction in nutrient loadings. Chlorophyll concentrations for the spring and summer months for the inshore boxes were adjusted using the multipliers for either the 25 or the 50% nutrient reductions (Fig. 5). Two simulations were compared: chlorophyll adjusted as in baseline and hypoxia included (present day (PD)) and chlorophyll adjusted for a 25% reduction in loadings and hypoxia included. A parallel set of simulations were performed for the 50% reduction in nutrient loadings.

Simulation Set 4: Offsetting Reduced Nutrient Loading with Reduced Hypoxia
The fourth set of simulations was built upon the third set (25 and 50% nutrient reductions with present day hypoxia) by determining how much would DO need to increase (hypoxia reduce) to offset the negative effects on croaker of less nutrients leading to lowered food. In a series of simulations, the hourly DO concentration was increased from present day values by adding 0.25, 0.5, 1.0, or 1.5 mg l −1 to each cell.
Increasing the DO concentrations resulted in a delayed onset of hypoxia in the spring and lower hypoxia areas during the end of July The model-generated hypoxic areas (maximum 10-day average and July), which was limited to four distinct values (i.e., the FVCOM years) under present day hypoxia conditions, that shifted to lower hypoxic areas as DO were increased (Fig. 6). For example, under present day hypoxic conditions for years 101 to 140 and using the maximum 10-day average, the model simulation had 8 years of normoxia (mild), 21 years at 15,831 km 2 (2010), 7 years at 18,399 (2001), and 4 years at 17,596 (2002). Using the late July average (rather than the 10-day maximum), the same simulation had 8 years of normoxia (mild), 21 years at 8889 km 2 , 4 years at 15,253, and 7 years at 14,670. Both sets of hypoxia areas (using 8, 21, 7, and 4) are shown plotted as present day in Fig. 6. These values coarsely covered the historical July values estimated from the annual summertime cruises (Obenour column in Fig. 6), except for no model values that exceeded 20,000 km 2 . The mode and range of hypoxic areas progressively decreased with increasing additions of DO (PD to 0.25 to 0.5 to 1.0 to 1.5 in Fig. 6). For example, the maximum 10-day area for the severe year of 2002 was 17,590 km 2 in present day, and then to 14,570 for +0.25, 11,790 for +0.5, 7210 for +1.0, and 3492 for +1.5 mg l −1 . These values appear as the frequency of 7 in Fig. 6 because that is how many times 2002 were used as severe year in years 101-140 of the simulations. The intermediate hypoxia year of 2010 occurred 21 times, and went from 15,830 km 2 in present day to 6260 km 2 with the largest DO addition of +1.5 mg l −1 . The mild year was normoxic in present day and was therefore unaffected by increasing DO concentrations; it repeatedly appeared as count of 8 at zero area. The overall effect of adding DO is seen by the average hypoxic area (over years 101-140) steadily decreasing with increasing DO additions (solid circles left to right in Fig. 6). This is a result of each of the 3 FVCOM years (with their respective frequencies of occurrence) having lower hypoxia areas as DO is increased.
This fourth set of simulations was designed to address the general question: How much would DO need to increase from present day (due to reduced nutrient loadings) to reach the Bbreak-even^level where the negative effects of reduced food were offset by the benefits of improved DO? We used two versions of this question that differed by the target selected to determine the break-even level. Both versions started with the reduced nutrient simulation (either 25 or 50%) with present day hypoxia. The first version used a target of croaker abundance simulated under present day conditions (from set 3), which was chlorophyll multipliers set to baseline (not reduced nutrients) and baseline or present day hypoxia (i.e., the first set of simulations). This target could be considered as the (albeit very unlikely) situation of croaker being the only consideration about the efficacy of reduced nutrients. How much does DO have to increase so that there is no loss, or even a gain, in croaker despite reduced nutrient loadings?
The second version used a target of croaker abundance under the reduced nutrient loadings but with normoxia every year (from simulation set 1). This target is the absolute highest abundance croaker could obtain under the reduced nutrient conditions (i.e., no hypoxia). Simulated croaker abundance can never exceed the values predicted under normoxia, so we are looking at what increase in DO results in croaker abundance approaching this target value. Use of reduced nutrients with normoxia as a target can be considered the situation when reducing nutrients are decided for other reasons than croaker, and approaching the target means that croaker will do as well as possible under the reduced nutrients.

Simulation Set 5: Sensitivity Analysis
The fifth set of simulations was a sensitivity analysis to explore the effects of low DO avoidance and the relative contributions of growth, mortality, and reproduction to the overall population response (Table 1). We used the two ways we represented hypoxia in the model: idealized (from Rose et al. 2017a) and from FVCOM (this paper). The only difference between these two sets of sensitivity simulations was the degree of variability of hypoxia (idealized versus FVCOM). We used the uncoupled version with FVCOM DO here because that allowed for a direct comparison to the results using idealized hypoxia reported in Rose et al. (2017a), which also used climatological chlorophyll.
The effect of avoidance was explored by reducing the number of cells searched in the neighborhood from 81 to 25 for all individuals and by using 1.5 and 2.5 mg l −1 DO values (rather than 2.0 mg l −1 ) to trigger the avoidance response. We also explored the cues used for selecting the destination cell as part of avoidance by ignoring temperature when selecting alternative cells, and eliminating the direct effects of low DO exposure on growth, survival, and reproduction (i.e., we set V G , V S , and V R to one) when temperature was ignored. The contributions of growth, mortality, and reproduction were explored by performing the hypoxia simulation with only mortality included and then with mortality and reproduction and mortality and growth. The contribution of reproduction or growth was then computed as the percent of Bmissing^abundance obtained when all three are included. To illustrate, suppose hypoxia simulation with all effects had 1000 fewer individuals in the 40-year average abundance than under normoxia, the simulation with mortality only resulted in 600 fewer individuals than with normoxia, and with mortality and reproduction had 800 fewer individuals. The contribution of reproduction would be estimated as [(800 − 600) / 1000] or 20%.
We included three more simulations in the sensitivity analysis to further explore the contribution of mortality, reproduction, and growth under FVCOM DO. Because the population response shifted from reproduction under idealized to direct mortality under FVCOM, we wanted to assess the robustness of predictions if we had overestimated the direct mortality effects. First, we reduced the no-effects concentration (Eq. 3) for the mortality effect (DO NE from 1.35 to 1.1 and to 0.8). Second, we reduced the DO NE for mortality to 0.8, but also increased the DO NE for reproduction (4.4 to 5.5) and growth (3.5 to 4.5). Changing the DO NE not only affects the DO concentration at which effects begin but also the sensitivity to all DO concentrations below DO NE because it changes the shape of the curve that determine the f values. Higher values of DO NE shifts the curve down, so higher multipliers (greater reductions) are obtained for any give DO below DO NE ( Fig. E-1, Rose et al. 2017a).

Model Outputs
The basic model output reported for all sets of simulations was the annual abundance of age 2 and older individuals for years 1 to 140. Age 2 and older abundance was computed on September 1 just after birthdays (i.e., youngest individuals just turned age 2). The annual abundances were further summarized by reporting the average abundance for model years 101 to 140 and also expressed as a percentage of the baseline also averaged for years 101 to 140.
Some additional model outputs were reported for simulation set 1 to help explain the reasons for the predicted responses of croaker to FVCOM-derived hypoxia. The outputs were as follows: mean stage survival of late juveniles, age 1, and age 2; averaged length at the end of the late juvenile (entering age 1), age 1 (entering age 2), and age 2 (entering age 3); mean eggs per batch per gram (EPG) and mean eggs per batch per individual (EPI) of age 1, age 2, and age 3; mean fraction of mature age 1 croaker; and total annual egg production. The details of how these are computed are provided in Table 5 of Rose et al. (2017a). Stage (and age) survival captures the direct mortality effects of low DO, and averaged length-at-age and fraction of age 1 that were mature would show growth effects. Eggs per batch per gram (EPG) illustrate the effects of reduced fecundity due to direct effects of low DO, while eggs per batch per individual (EPI) show the direct effect of DO on fecundity plus any effects due to slowed growth and associated smaller size at age. Total annual egg production reflects the combined effects of growth, survival, and reproduction. We computed averaged values of each output variable for all hours in each year of the simulation, and grouped the year-specific values according to whether they were a mild, or 2001mild, or , 2002mild, or , or 2010 To isolate the effects of hypoxia, we compared the hypoxia versus normoxia for uncoupled chlorophyll and coupled chlorophyll; the only difference between the pairings was hypoxia or not. Averaged values of the outputs were expressed as the ratio of the average values in hypoxia to normoxia.

Isolated Effect of Hypoxia (Set 1)
Inclusion of hypoxia resulted in about a 25% reduction in long-term croaker abundance for both chlorophyll uncoupled and coupled to hypoxia (Fig. 7). Coupling chlorophyll to Table 1 Summary of sensitivity runs of the croaker IBM showing the factor examined, the labels used in Fig. 11, and details of how the simulations were performed Label used in Fig. 11 Factor examined Details severe year and, for survival, had larger effects than the 2002 severe year. Survival of the late juvenile stage increased by 7 to 10% because lower abundances lead to reduced mortality via the density-dependent effect. Fraction surviving age 1 and age 2 were the same as baseline for mild years (as expected) and the 2002 severe years, and about 5% lower for the intermediate and 12% lower for the 2001 severe years. These combined to result in a 5% lower survival when averaged over the entire 40 years. Average length entering age 2 and age 3 and the fraction of age 1 mature (length dependent) were unchanged from baseline (normoxia) for all years and for coupled and uncoupled chlorophyll. EPG and EPI showed only slight reductions for age 1 individuals. Age 1 had minimal exposure because they spawned in their first year and thus were exposed mostly as late juveniles that were restricted to depths less than 20 m. EPG and EPI of age 2 and age 3 showed similar reductions as survival except for the 2002 severe years. EPG and EPI were not reduced in mild years, reduced by 6-7% in the intermediate years, and by 10-11% reductions in the 2001 severe years. However, unlike survival, the 2002 severe years showed similar reductions in EPG and EPI as the 2001 severe years (10-13%), rather than the 2002 severe years showing no effect as with survival. Total annual egg production (lagged 1 year forward to account for summer exposure affecting spawning next year) was reduced by about 28% for all years and simulations with hypoxia. The consistent effect across mild to severe years is because egg production integrated effects across age classes that had experienced a different mix of hypoxia year types.

Mortality of Benthos (Set 2)
There was no croaker population response under the realistic assumption of a direct mortality effect on croaker food (Fig. 8).
The lack of the effect of reduced food was due to only marginal reductions predicted in the P-values and because these reductions occurred in cells with DO less than 2 mg l −1 which also were effectively avoided by croaker. When the mortality effect was unrealistically increased by reducing the value of Q (i.e., from 1000 to 100 h), the model did predict a substantial population response showing population responses are possible with our representation but just not realistic.

Reduced Nutrient Loadings and Offsetting (Sets 3 and 4)
The 25 and 50% reduced nutrient loadings (with present day hypoxia) resulted in further decreases in croaker abundance due to lower food (blue versus red, Figs. 9 and 10). Both of these simulations included the same present day hypoxia conditions, but differed only in their nutrient loadings (i.e., chlorophyll multipliers). The 25% reduction in nutrient loadings resulted in an average abundance of age 2 and older of 6.32 million compared to 8.08 million under present day conditions, a reduction of 22% due to lowered food. The 50% reduction in nutrient loadings resulted in a 47% reduction in croaker abundance from present day food conditions. For the 25% reduction in nutrient loadings, an increased DO concentration of 0.5 to 1.0 mg l −1 would result in no loss to croaker (Fig. 9). The target of present day abundance was a lower abundance (easier to reach) than the abundance with a 25% nutrient reduction with normoxia because of the relatively large positive effect of reducing hypoxia compared to the smaller negative effect of a 25% reduction in nutrients (i.e., normoxia (green) above present day (blue) in Fig. 9). The present day abundance was between the simulations under a 25% nutrient reduction with DO increased by 0.5 and 1.0 mg l −1 . Thus, croaker would show no loss by the 25% reduction in nutrient loadings if the reduced nutrients were associated with a 0.5 to 1.0 mg l −1 increase in DO. A 1.0 mg l −1 increase also resulted in croaker abundance approaching the other abundance target of a 25% reduction with normoxia (other reasons than croaker for reducing nutrients). Such a decrease in hypoxia under a 25% reduction in loadings would result in croaker abundance being as high as is possible with the lowered food.
For the 50% reduction in nutrient loadings, an increased DO of 1.0 mg l −1 would result in croaker doing as well as possible +0.25 Fig. 9 Annual age 2 and older croaker abundance for normoxia with reduced food (Normoxia, green), FVCOM-derived present-day hypoxia with baseline food (PD, blue), present-day hypoxia with reduced food (25% Nutrients, red), and DO adjusted by adding +0.25, 0.5, 1.0, and 1.5 mg l −1 to present day hypoxia with reduced food. Reduced food was due to the 25% reduction in nutrient loadings. Note the y-axis starts at 5 × 10 6  Fig. 10 Annual age 2 and older croaker abundance for normoxia with reduced food (Normoxia, green), FVCOM-derived present-day hypoxia with baseline food (PD, blue), present-day hypoxia with reduced food (50% Nutrients, red), and DO adjusted by adding 0.5, 1.0, and 1.5 mg l −1 to present day hypoxia with reduced food. Reduced food was due to the 50% reduction in nutrient loadings. Note the y-axis starts at 3 × 10 6 (normoxia target), but DO cannot increase enough for croaker abundance to approach their abundance under present day food conditions (Fig. 10). This is seen in Fig. 10 as the +1.0 mg l −1 line approaching the green line (normoxic with reduced nutrients) but remaining below the present day line (blue). The reduction in food under 50% nutrient loadings caused a larger negative effect than can be offset by the positive effect of eliminating hypoxia. Thus, if a 50% reduction in nutrient loadings resulted in a 1.0 mg l −1 increase in DO, croaker was predicted to do as well as possible with the reduced food (reductions in loadings justified for other reasons) but, with croaker viewed as the sole reason for reducing loadings, there would be no benefit to croaker. The simulation with reduced nutrients and improved hypoxia produced less croaker less than present day conditions.

Sensitivity Analysis (Set 5)
Model predictions of age 2 and older croaker abundance were more sensitive to the size of the neighborhood searched in avoidance (A in Fig. 11) than to the DO concentration that triggered avoidance (B and C), with the FVCOM derived (dynamic DO) causing larger reductions than the idealized hypoxia (top versus bottom panels in Fig. 11). Reducing the cells in the search neighborhood from 81 to 25 (A) resulted in a hypoxia effect of 0.8 (versus 0.89 for present day, shown as line) for idealized and 0.41 (versus 0.72) for FVCOM DO. Under idealized hypoxia, the smaller neighborhood delayed how long individuals needed to find higher DO concentrations but the spatially smooth and temporally stable DO concentrations with idealized hypoxia allowed individuals to still effectively avoid low DO. With DO varying under FVCOM hypoxia, the smaller neighborhood made it even more difficult for individuals to avoid low DO for any sustained time.
Only changing the threshold DO that triggered avoidance (1.5 or 2.5 versus 2.0 mg l −1 in present day), resulted in relatively small changes in average age 2 and older abundances (B and C in Fig. 11). Predicted population reductions for 1.5 (B) and 2.5 (C) mg l −1 triggers were 0.73 and 0.89 versus 0.89 for present day for idealized (top panel in Fig. 11) and 0.64 and 0.68 versus 0.72 for FVCOM DO (bottom panel in Fig. 11). The larger effect of the 1.5 mg l −1 threshold under idealized hypoxia was due to allowing individuals to be exposed to DO concentrations between 1.5 and 2 mg l −1 , which had significant effects on growth and reproduction of individuals. The larger effect of a reduced neighborhood size under FVCOM DO carried over when combined with the different thresholds. Any simulations with different triggers that included the smaller neighborhood with FVCOM DO resulted in relatively large reductions in population abundance (D and E in Fig. 11).
While variability in DO resulted in higher sensitivity to neighborhood size and threshold values, the effect of selecting temperature for the destination cell was more important under idealized hypoxia and therefore indirect effects on growth played a more important role under idealized conditions. When we changed from individuals selecting cells within their neighborhood using the optimal temperature to individuals selecting cells with highest DO (ignored temperature), the effect of hypoxia was a lower average population under idealized (0.75 versus 0.89) but a similar reduction in abundance (0.72 versus 0.75) under FVCOM (F in Fig. 11). The varying DO under FVCOM hypoxia eliminated any costs of going to less optimal temperatures because the dynamic DO prevented any locations from having both sustained tolerable DO and near-optimal temperatures. Using DO as the destination but eliminating direct effects of low DO (V S = V R = V G = 1, G in Fig. 11 Table 1 with direct effects). Thus, when DO conditions are smooth in time and space, the cue used to select the destination cell under avoidance movement is important and the indirect effects (temperature at destination) of displacement can play a relatively important role. With varying DO concentrations, direct effects of low DO dominated the larger population response. Reduced reproduction dominated the population response under idealized conditions, whereas increased mortality, and to a lesser extent reduced reproduction, combined to determine the population response under FVCOM hypoxia (H, I, and J in Fig. 11). The consistent DO conditions under idealized conditions enabled individuals to avoid lethal conditions but allowed for sublethal exposures with reproduction that was more sensitive than growth. For example, at a DO concentration of 3.0 mg l −1 , the f R value for reproduction was 0.63, while it was nearly no effect (0.97) for growth. Under FVCOM hypoxia, survival alone (H) generated a majority of the population reduction (0.8 versus 0.72), and when combined with reproduction (J) generated the entire response (0.75 versus 0.72). Using the mortality, mortality with reproduction, and mortality with growth simulations (H, I, and J in Fig. 11), approximately 62% of the response under idealized hypoxia was due to reproduction effects, 26% due to slowed growth, and 12% due to direct mortality. These shifted under FVCOM DO, with mortality accounting for about 70%, reproduction about 26%, and growth only about 8%.
The sensitivity analysis simulations with FVCOM DO only that varied the no-effects parameters (DO NE ) to explore the degree of dependence of model results on direct mortality effects showed that the population responses to hypoxia were robust to the heavy reliance on mortality effects (labeled Changed Sensitivity in Fig. 11). When we just lowered the DO NE for mortality (1.35 to 1.1 mg l −1 , K in Fig. 11), the effects of reproduction and growth increased so that the same population reduction (0.76) was predicted. Even when the DO NE for mortality to 0.8 mg l −1 (L), which resulted in very little direct mortality, the effects of reproduction and growth, now applied to individuals who had previously died under with higher DO NE value for mortality before, still resulted in a 15% reduction in population abundance. Increasing the DO NE for growth (4.5 from 3.5) and reproduction (5.5 from 4.4), coupled with a lowered DO NE for direct mortality, resulted in a similarly large population reduction (0.68 versus 0.72) as with the original DO NE values (M in Fig. 11). In this case, the contributions were shifted to 46% from reproduction, 42% from growth, and only 19% from mortality. Survival of age 1and age 2 adults (averaged for all years 101-140) showed a smaller reduction compared to the original present day (PD) hypoxia results (0.98 versus 0.95, Table 2), while EPI of age 2 and age 3 showed larger reductions (0.82 versus 0.94); weight-at-age, which showed no effects with the original DO NE values, now showed a 4% reduction for age 2 and age 3.

Discussion
Our best estimate based on our modeling results is that hypoxia in the northern Gulf of Mexico reduces long-term average croaker population abundance by a minimum of 10% and most likely by about 25%. Rose et al. (2017a) presented a minimum estimate of a 10% reduction that used idealized hypoxia conditions and equipped the simulated croaker individuals with very effective avoidance behavior. There was some direct mortality during the 1-week setup period of the idealized hypoxia, but then the individuals avoided lethal concentrations after that because the restricted area search was very effective when low DO concentrations remained relatively constant over time and space. In this paper, we used more realistic DO conditions generated by a 3-D hydrodynamicswater quality model, included the situation of croaker food being correlated to the severity of hypoxia, and performed a sensitivity analysis on the avoidance movement algorithm and other assumptions.
When we used DO concentrations from the FVCOM model for three different years strung together into a 140-year time series, we predicted about a 25% decrease in long-term (final 40 years) croaker abundance. Coupling croaker food with hypoxia resulted in not only higher population abundances but also still about a 25% reduction due to hypoxia. Sensitivity analysis that used DO from FVCOM showed that the variation in avoidance assumptions often also generated reductions of 25% or more. The detailed analysis of the hypoxia simulation using FVCOM output and the sensitivity analysis also showed that mortality became more important than under the idealized hypoxia conditions because DO varied hourly thereby creating very low but persistent episodic exposures to lethal concentrations. Effects on fecundity remained important under both idealized and FVCOM (dynamic) DO conditions because exposure occurred above the avoidance trigger so variation in DO did not greatly affect the already relatively high exposure to sublethal DO under idealized conditions that results in fecundity effects. There is substantial field evidence that hypoxia has a negative effect on the reproductive output of croaker in the NWGOM (Thomas et al. 2015;Rahman 2009, 2012). Recent measurements also showed a potential reproductive effect of hypoxia on croaker in the Chesapeake Bay (Tuckey and Fabrizio 2016). We use a 25% reduction in abundance (rather than higher) as our most likely estimate because our use of only of a few years from FVCOM limited how well we could represent intermediate and severe hypoxia years and the larger role played by direct mortality that can only be crudely represented on the kilometer and hourly resolution of the IBM.
Subsequent simulations show that our estimate of a 25% reduction was robust to the role played by direct mortality. The effects of mortality, reproduction, and growth are dependent on each other in the model; those individuals that do not die from low DO will go on and exhibit reduced fecundity and growth. When we lowered the sensitivity to mortality by lowering the value of DO NE , similar population reductions were still predicted because those that survived showed even higher reproduction and growth effects. Further, when we decreased sensitivity to mortality and increased sensitivity to growth and to reproduction (i.e., shifted to 80% of the population response due to reproduction and growth), we predicted similar or larger than 25% population reductions. The changes in DO NE values we used (Table 1) were within the realm of possibility given the variability in the results of the specific laboratory experiments we used (Fig. ESM-2 in Rose et al. 2017a) and additional variability typically observed across experiments (e.g., EPA 1986EPA , 2000 and in aspects of field avoidance (e.g., Wannamaker and Rice 2000;Craig 2012;Craig and Bosman 2013). While we are less certain of the mix of mortality, reproduction, and growth effects, we have confidence that the combined effects of many possible realistic mixes support our determination that our best estimate is about a 25% population reduction.
We focused on low DO directly affecting growth, mortality, and reproduction of individuals and on low DO reducing prey of croaker. The direct effects on mortality and growth are well documented in laboratory experiments; the effects on reproduction are not typically considered when discussing hypoxia effects on fish and they played a role in the population responses simulated here. The simulation of direct mortality from low DO reducing croaker prey had no effect on croaker population abundance. Croaker appears to remain benthic feeders under hypoxia rather than switching to a more pelagic diet (Mohan and Walther 2016). However, there is some evidence that, at least under episodic hypoxia, croaker can even, at least temporarily, increase their predation rates by taking advantage of more vulnerable prey (Long and Seitz 2008). There are other pathways by which hypoxia could affect croaker population dynamics that were not included in our model simulations. For example, hypoxia avoidance can lead to changes in overlap of croaker and their predators (see Ludsin et al. 2009), hypoxia exposure can reduce sperm production by males (Thomas and Rahman 2010), and shifts in shrimp fishery in response to hypoxia (Smith et al. 2014) can affect croaker mortality via bycatch. These additional effects can be added to the model simulations (e.g., predation from habitat compression- Adamack et al. 2017). We focused on direct effects on growth, mortality, and reproduction (and reduced prey) because they had data available, and laboratory and field results showed them to be important at the individual level and therefore potentially important at the population-level.

Evaluation of Several Critical Assumptions
There are three assumptions that are both critical to the results and can be considered highly uncertain. These are the following: no vertical avoidance movement, assignment of DO values to the relatively large 1-km 2 cells every hour, and use of surface chlorophyll concentration as an index of croaker food. We present evidence for using a 2-D horizontal approach based on field data and additional avoidance modeling. We also provide modeling results and field data on croaker exposure to hypoxia to show that not allowing vertical movement and using 1-km 2 cells, even if only crude approximations, resulted in realistic exposure in the model simulations. Finally, justification for using surface chlorophyll as an indicator for croaker food is based on field studies showing that, to various degrees, primary production influences benthos abundance.
Evidence supporting a 2-D horizontal approach for croaker includes both field data and additional avoidance movement modeling. Craig (2012) reported the results of comparing trawl catch rates during summer of 2002 and 2004 between tandem nets such that the lower net remained in hypoxia water while the upper included the hypoxic layer but also extended just above the hypoxia layer. The height of the lower net averaged 1.3 m off the bottom and was in <2 mg l −1 water for 84% of the time. In contrast, the upper net was located above the 2 mg l −1 by an average distance of 2.2 m for 80% of the time. Craig (2012) reported that there were no significant differences in croaker catch rates between the two nets regardless of DO concentrations, suggesting that croaker were not moving up vertically when DO was low. He did find differences between the two nets related to DO in other species confirming that the data and analysis could detect such vertical avoidance. Horizontal movement (displacement) of croaker in response to hypoxia has also been documented from field data (Craig and Crowder 2005;Craig 2012).
LaBone (2016) used the same FVCOM-WASP model as used here and, in one analysis, focused on the effects of allowing vertical avoidance behavior on croaker exposure to low DO concentrations. She used a restricted area search algorithm for horizontal and vertical avoidance movement, as we used here in the IBM for horizontal avoidance only, and two random-walk-based algorithms for other behaviors. Simulations using about 1000 individuals at 15 min time steps were done for July 24 through August 2, 2002, with all individuals confined to the bottom layer (2-D horizontal movement only) and alternatively also allowed to move vertically in response to low DO (i.e., 3-D). When the individuals were error free in their movement decisions, all combinations of avoidance and default movement algorithms that included vertical movement resulted in much reduced exposure to very low (<2 mg l −1 ) DO. With increasing rates of error or misinformed movement decisions, exposure to very low DO increased with vertical movement but was always substantially less than for the corresponding 2-D simulation. However, exposure to DO concentrations of 2 to 4 mg l −1 (growth and fecundity effects in our IBM model) was similar, or even at times higher, when vertical movement was allowed. These results suggest that allowing croaker to move vertically and relatively error free would reduce the direct effects related to mortality but would not reduce the direct effects due to fecundity and growth. In the IBM, it was individuals exposed to nonlethal DO concentrations that accounted for most of the reductions in growth and fecundity because increased mortality effects occurred quickly and thus those individuals exposed to lethal concentrations died rapidly and thus played a small role in growth and fecundity effects.
The combination of no vertical movement and a relatively coarse grid (1-km 2 cells at 1-h time steps) in the IBM could result in individuals being over-exposed to low DO concentrations. The minimum exposure is necessarily 1 h and avoidance behavior must be effective enough to move individuals kilometers to get to new cells and higher DO concentrations. Simulated exposure in the IBM was realistic. We analyzed the IBM simulation with FVCOM DO and chlorophyll multipliers under present day conditions ( Fig. 9) to determine exposure levels. We focused on the two severe hypoxia years (2001 and 2002) because this was when over-exposure was most likely to occur. Using the same approach as in Rose et al. (2017a), we computed the percent of all age 0 to age 2 (i.e., all vulnerable) individuals grid-wide that were exposed to <1.35, 1.35-3.35, and 3.35-4.4 mg l −1 for each hour in year 138 (2001) and year 136 (2002). While the percent exposed varied from day to day as hypoxia changed, typical maximum hourly percent of all age 0 to age 2 exposed to <1.35 (lethal) was infrequent and <1% when it occurred. Maximum values of exposure were 10-25% for 1.35-3.35 and for 3.35-4.4 mg l −1 . Further, about 50% of the age 0 to age 2 individuals experienced at least 1 h of DO less than 4.0 mg/L during the summer in both 2001 and 2002 (i.e., 50% of the individuals never experienced any low DO).
We re-analyzed the fecundity data reported in Thomas et al. (2015) to provide a basis for evaluating the 10-25% exposure levels predicted in the IBM. Croaker was sampled over a broad area of the NWGOM in the fall that covered the entire area of known hypoxia and was at stations on transects (C, D, F, H, J, and S) used for the annual summertime hypoxia cruises (see Rabalais et al. 2001). In addition, croaker was sampled in the eastern part in the Mississippi sound as a reference site. Here, we divided each individual's fecundity by the average fecundity from the reference site. Fecundity at the reference site could differ from fecundity in hypoxiainfluenced sites by more than just hypoxia. Using the reference site fecundity as normoxia fecundity is reasonable because the reference site fecundities were very similar to fecundities measured at normoxic sites and reported for Gulf of Mexico croaker (Thomas et al. 2015). The cumulative fraction of individuals was plotted against the normalized fecundity (Fig. 12).
To provide some idea of what kind of exposure would result in such reductions in fecundity, we used the exposureeffects model (based on laboratory exposures of 10 weeks) as used in the IBM. The vertical lines in Fig. 12 show how different combinations of exposure in the laboratory within a 10-week period would reduce fecundity. Figure 12 shows that the vast majority of individuals spread over a large area of the NWGOM in areas of hypoxia were exposed to the lab equivalent of weeks of low DO. For example, over 75% of the individuals were exposed to the laboratory-equivalent exposure of at least 2 mg l −1 for 10 weeks. Over 90% of the individuals showed exposure that is the equivalent of 3 mg l −1 for 5 of the 10 weeks (or 3.5 mg l −1 for 10 weeks) or worse. Similar results and high levels of exposure (most all individuals and weeks of low DO) were also measured in 2006 and 2007 at fewer transects. In addition, hypoxia-inducible-factor-1 α (HIF-1α) mRNA level was measured in 2007 and 2008 in croaker in the summer months (Rahman and Thomas, in review) and, based on extensive laboratory experiments with known exposures, showed most all individuals in the field had been exposed to low DO (<~3.5 mg l −1 ) for days including within the recent few days (see review in Thomas and Rahman 2009). This field evidence shows that the simulated exposure in the IBM, despite no vertical movement, hourly time steps, and 1-km 2 cells, was reasonable.
The use of surface chlorophyll concentration as an index of croaker food availability is only a first order approximation but has some empirical support. A preferred approach would The vertical lines show the fraction of normoxic fecundity that would results in 10-week lab exposures that were used to derive the exposure-effects submodel used for hypoxia effects on fecundity in the IBM (see Rose et al. 2017a) be to use the output of a 3-D model that simulated water quality (temperature, salinity, DO) as well as zooplankton and benthos dynamics for the NWGOM; however, such a model is not presently available as the modeling focus has been on hypoxia itself or on the food web responses but not an integrated model of water quality, the plankton food web, and its coupling to the benthos (Rose et al. 2017c). We therefore used chlorophyll as a rough proxy for croaker food.
In general, across system comparisons show that chlorophyll, or primary production, is clearly related to fish productivity, including pelagic and demersal fish species, and the most likely explanation is increased food resources (Nixon and Buckley 2002;Caddy 2000;Breitburg et al. 2009). Where analyses have been conducted, positive relationships between chlorophyll and fish abundance or landings (including demersal species) are reported across subregions within coastal systems (not just across systems) as well (Blanchard et al. 2012;Ware and Thomson 2005;Frank et al. 2006). What is needed for the IBM is an indicator that captures the relative variation in food available to croaker, not the absolute magnitude. We calibrated chlorophyll concentrations to the P-value from the bioenergetics model (Fig. 2), thereby removing the dependence on absolute values.
An assumption of our approach is that chlorophyll is an indicator of food for croaker on a regional and seasonal basis within the NWGOM. While juvenile and adult croaker are clearly benthic feeders, consuming a variety of benthic meiofauna and macrofauna (Overstreet and Heard 1978;Nye et al. 2011), large-scale benthos sampling, with comparable measurements of primary production, is limited in the Gulf of Mexico (see data and review in Rowe and Kennicutt 2009). The Mississippi River seasonally enhances nearshore surface primary productivity Dagg and Breed 2003). Spatial patterns in the flux of particulate organic carbon to the sediment parallel patterns in surface productivity, indicating a coupling between riverborne nutrients, surface phytoplankton, and organic carbon inputs to the sea floor (Wysocki et al. 2006;Green et al. 2008;Rabalais et al. 1996). This elevated vertical flux of organic matter likely enhances benthic productivity in riverine-influenced regions of the shelf, though direct measurements are lacking. Surface chlorophyll-a concentrations have been correlated with the summer (June-July) spatial distribution of cownose rays (also a benthic feeder like croaker) at spatial scales of kilometers across the northwestern GOM shelf (Craig et al. 2010).
More debatable is whether this relationship between chlorophyll and benthos productivity persists further along the coast away from the plume influence and into deeper waters on the shelf. Baguley et al. (2006) used biweekly average values for November 1999 through April 2000 and showed that meiofauna abundance increased with increasing chlorophyll-a concentration in the water column. They stated BDespite considerable heterogeneity in primary production in northern Gulf of Mexico surface waters; however, a relationship with benthic meiofauna abundance exists.^The chlorophyll used by Baguley et al. was from SeaWiFS, which we also used to estimate the chlorophyll values in the IBM (Rose et al. 2017a). A similar positive correlation between water column chlorophyll and benthos (25 stations spread over 4 transects) was reported off of Texas (10 to 134 m) on seasonal transects extending from shore to the shelf for 1976-1977(Flint 1981. Thus, use of chlorophyll as an indicator of food availability appears valid for the NWGOM as long as it is used on sufficiently coarse spatial and temporal scales. The connection between chlorophyll and benthos is neither direct nor strong enough to allow for chlorophyll to be used to resolve sub-kilometer and weekly scale variation in food availability. Given the limited data, our use of interpolated (smoothed) chlorophyll as an indicator of food availability on a monthly basis for our grid (Rose et al. 2017a) and then for spring and summer seasons in the nearshore region of the Gulf under nutrient reduction loadings scenarios appears appropriate.

Feasibility of Break-even Points
We accumulated predictions of hypoxia area under reduced nutrient loadings to compare to the area reductions projected for the IBM under increasing DO concentrations (Fig. 13).
The purpose was to determine if the identified additions of DO (and resulting reductions in hypoxia area) that offset the croaker losses due to less food were possible. For example, we identified that a DO increment of +0.5 mg l −1 was sufficient to offset croaker losses due to less food under a 25% reduction in nutrient loadings. The question address here is whether the reduced hypoxia area under a 0.5 mg l −1 increase in DO was feasible with a 25% reduction in nutrient loadings.
We used the published statistical model results for predictions of late July area for reduced nutrient loadings. These models all used the historical time series of July hypoxia areas (either original LUMCON estimates or re-estimated in Obenour et al. 2013) and nutrient (nitrogen or phosphorus) loadings for the Mississippi and Atchafalaya Rivers. The statistical models differed in the percent reductions of the nutrients simulated, and also in their explanatory variables (e.g., nitrogen alone, nitrogen and phosphorus) and formulations (e.g., transformations, lagged terms). One set of model predictions we used  were reported as 5year averages of predicted areas. We reported each of the statistical model predictions as the fraction of hypoxia area for each historical year (or 5-year average), with the predicted area under nutrient reduction divided by the historical area (no nutrient reduction) for the same year (or same 5 years). We used year-specific ratios as reported by each statistical analysis to show the spread in predicted area reductions and because the statistical models and the IBM did not all start with the same historical values of areas.
The break-even values of increased DO, and corresponding July hypoxia areas, identified for the 25 and 50% nutrient reductions appear to be feasible. The increase in DO of 0.5 mg l −1 (needed for the 25% nutrient reduction approaching present day) was centered on about 0.75 (solid horizontal line in Fig. 13), which corresponded to an average of a 25% reduction in the hypoxic area. All of the statistical models and FVCOM, except one model, predicted ratios that were lower or surrounded the reduced area values for the 0.5 mg l −1 DO increase in the IBM. This is seen in Fig. 13 as the solid line (+0.5 mg l −1 ) being within or above most of the points on the left-side panel. The one exception was Obenour et al. (2012) without sediment oxygen demand (SOD) included (average of 0.87 or a 13% reduction). Feist et al. (2016) used a water quality model (not included in Fig. 13) and also reported about a 22% reduction in hypoxia area over 5-years with a 25% reduction in N and P loadings; however, Scavia et al. (2013) predicted that a larger roughly 40% reduction in nutrients was needed.
The likelihood of a 25% reduction in nutrients resulting in a 1.0 mg l −1 increase in DO (exceed present day and approach normoxia) is low. An increase by 1.0 mg l −1 resulted in an about a 50% reduction in hypoxic area (0.47) in the IBM (dashed line in left-side pane of Fig. 13). Three of the statistical model predictions under a 25% reduction in nutrients were less or similar to the 0.5, but these were predictions from the same statistical model  with different nutrient reductions. The remaining statistical models and FVCOM predicted that a 25% reduction in nutrients reduced areas would not achieve a 50% reduction in hypoxia area (most all values were greater than 0.5, or less than a 50% reduction).
The break-even value of increased DO of 1.0 mg l −1 for the 50% reduction in nutrients (no loss-normoxia target) was well within the values predicted by the statistical models and FVCOM but not ROMS (right-side panel in Fig. 13). Whereas it was very unlikely that a 25% reduction in nutrients could result in a 50% reduction in hypoxia area, a 50% reduction in nutrients is more likely to be able to generate such a reduction in hypoxic area. Three different statistical models included predicted areas that averaged less than 0.5 (i.e., more than a 50% reduction); the only exception was again the Obenour et al. (2012) prediction without SOD included (average of 0.75 or only as 25% reduction). Most points were below the +1.0 mg l −1 line in the right-side panel of Fig. 13, suggesting at least a 50% reduction in area was likely. While FVCOM predicted 0.54, similar to 0.47 estimated by adding 1.0 mg l −1 DO within the IBM, and which agrees with a simpler but spatially explicit model by Scavia et al. (2013), ROMS predicted a much smaller reduction in area for 2004 to only 0.8 (20% reduction in area).
Converting the reductions in hypoxic areas to absolute areas so the IBM results of croaker abundances can be related to hypoxia areas is tenuous without further simulation analyses that use many more FVCOM years to replay the range of conditions observed in the historical record. At this point, we can conclude that (1) reductions in hypoxia needed with a 25% reduction in nutrients to result in no loss (or only moderate losses) of croaker is feasible and (2) croaker will do as well is possible (approach normoxic abundance) under the 50% reduction in nutrients.

Modeling Strategy
Our strategy for modeling was to use a bottom-up approach to ensure the key components of the coupled models were credible and realistic, and to then combine the components in order to simulate population-level responses to hypoxia. This involved modeling, laboratory experiments, and field data collection . This is a different approach than the usual fish population modeling approach whereby data on population abundance indices are used to calibrate and then validate the population model by comparing predicted versus observed time series (Jackson et al. 2000;Stow et al. 2009). We opted for a bottom-up approach to modeling because many versions of the model and a wide range of realistic parameter values could generate output that would be within the variability of empirical indices related to croaker abundances. Likely, versions of the model could be found that include hypoxia effects in various ways that would fit the data equally well. We therefore carefully developed the components of the modeling using existing and purposely designed experiments and field data collection to ensure credibility of predicted population responses. While there are many descriptions about how hypoxia could potentially affect fish at the population and higher (community, food web) levels, there are fewer quantitative estimates of these effects. One group of studies uses conceptual and logical arguments to infer that population-level effects are feasible and should be likely in some cases. A diagram, originally proposed by Caddy (1993), that shows increasing fish productivity or landings with increasing nutrient loadings until hypoxia effects become important and productivity or landings then decrease, has appeared in many papers (e.g., Diaz and Solow 1999;Caddy 2000;Diaz 2001;Hondorp et al. 2010). These papers also often go on to cite a few case studies, such as the Baltic or Black Seas, as examples where fisheries decreased due to eutrophication and hypoxia (e.g., Caddy 2000;Diaz 2001;Breitburg et al. 2009;Ekau et al. 2010). Many of these reviews conclude that there are local effects of hypoxia on fish (sometimes involving many individuals), but there is a general lack of evidence for hypoxia effects on fish populations (Wu 2002;Wannamaker and Rice 2000;Breitburg 2002). Breitburg (2002) cites several examples where hypoxia affected recruitment, but then concludes that BThe brevity of the list of clear examples of large-scale, hypoxia-induced population declines likely also reflects the conflicting effects of decreased habitat extent and suitability, and increased prey production….Â second type of analysis compares empirical data on fish abundance or landings and nutrient loadings across ecosystems to infer how a single system may respond over time to eutrophication (space for time). Caddy (2000) plotted nutrient availability versus the ratio of landings of pelagic to demersal species (P/D) for 11 semi-enclose marine areas. He noted a decrease in demersal relative to pelagic fish landings, although only as far as nutrient loadings were inferred to occur with seasonal hypoxia; there were no values in the region of the plot where year-round hypoxia (far right-very high nutrient loadings) was inferred to occur. A more recent update of this approach also attempted to determine the correlation of fish abundances with hypoxia. Breitburg (2002) noted a less than expected increase in demersal fish production (relative to pelagic fish prediction) for 14 semi-enclosed European Seas, which suggested a negative effect of hypoxia; however, there was only a weak correlation between harvests of demersal species and several indicators of hypoxia severity. Hondorp et al. (2010) analyzed 22 ecosystems and found a relationship with increasing P/D with hypoxia, but a negative relationship of landings of demersal species with hypoxia depended on the inclusion of single point (Black Sea). They discuss the importance of site-specific fishery management practices on affecting variability in fishery landings. Breitburg et al. (2009) analyzed 30 systems that overlapped with the analysis of Hondorp et al. (2010), and concluded that there was not much support for negative effects of hypoxia that were separable from the effects of increased nutrient loadings. A different approach was used by Chesney and Baltz (2001), who analyzed the time series of landings and catch-per-unit-effort (CPUE) for the NWGOM and concluded that Bthe exploited nekton populations in the Gulf of Mexico are able to tolerate the effects of hypoxia without obvious major consequences to their recruitment, production or population health.^In a third paper (Rose et al. 2017b), we used the 25% reduction in croaker abundance and used field data to add interannual variability to the predicted recruitment. We then tried to detect statistically the known reduction and demonstrated that detecting the effect (knowing that is was present) would be difficult without extensive sampling or information on recruitment dynamics. These types of correlative studies not finding obvious hypoxia effects have weak statistical power and thus leave open the question of whether hypoxia has population-level effects on coastal fish populations.
A third approach, and the approach most similar to the approach used here, is to focus on a single well-studied system (i.e., case studies) that has empirical data and to use modeling to isolate the hypoxia effect. Rose et al. (2009) attempted to synthesize the results from a suite of models that all focused on quantifying the effects of hypoxia from a single individual to populations within simple food webs. They concluded that evidence for large widespread direct effects of hypoxia was limited but that there were suggestions of situations in which direct and indirect effects on individuals could lead to effects on higher levels (e.g., cohort, recruitment, population). Eby et al. (2005) used measured growth rate reductions due to hypoxia observed in the laboratory and field, and then used a density-independent stage-based matrix projection model for croaker to estimate that hypoxia caused about a 4% reduction in population growth rate. This reduction, when extrapolated to 5 consecutive years of hypoxia, would translate into a 15% reduction in population abundance. Zhang et al. (2014) used a habitat approach with potential growth rate that depended on temperature, prey availability, and DO and concluded that the effects of hypoxia on total habitat quality in the northern Gulf of Mexico centered on Louisiana were small for bay anchovy and menhaden. Using EwE model for the Baltic Sea, Hansson et al. (2007) showed that hypoxia effects on reproductive volume of cod, if essentially eliminated, would result in about a 100% increase in cod catches over present day low values and catches similar to historically higher catches. Teschner et al. (2010) added an hypoxia effect to cod feeding and therefore cod predation on cod (cannibalism), sprat, and herring within a multispecies stock assessment model (virtual population analysis). They demonstrated that while cod biomass showed little difference without and with the hypoxia effect included, annual sprat and herring biomasses increased and decreased by 20% (at times up to 100%) in various time periods in the subareas analyzed. Finally, de Mutsert et al. (2016) recently used EcoSpace and concluded that under increased nutrient loadings, total fish biomass benefitted more from the increased primary production than the negative effects of more hypoxia, with the caveats that responses were specific-specific with a few species showing net loss of biomass. Model evidence for population or higher level effects of hypoxia on coastal fish species suggests that there are not widespread effects but that moderate to large effects are likely present for some species in some years.
While our final prediction of a 25% reduction in the longterm croaker abundance may not be surprising, it goes beyond circumstantial or conjecture and is defensible because of our methodology. Our approach was designed to isolate the population level effects of hypoxia on a well-studied species known to be exposed to sublethal effects of DO. The model development effort spanned multiple years in order to initially identify critical data and modeling gaps, and then fill these gaps with laboratory and field data collection and development of credible submodels (avoidance, exposure effects). An example of the data collection was the performance of laboratory experiments allowing for cause and effect relationships between low DO and croaker fecundity (e.g., Thomas et al. 2006;Rahman and Thomas 2007), and field studies documenting that the effects on fecundity observed in individuals in the laboratory was occurring at levels that could affect the population (e.g., Thomas et al. 2007Thomas et al. , 2015; Rahman 2009, 2012). A parallel effort in modeling was the development of an exposure-effects submodel designed for dealing with fluctuating exposures and that could be calibrated using widely available constant exposure lab experiments and then tested using the relatively few experiments that used fluctuating exposures (Neilan and Rose 2014). The robustness and credibility of our model-based predictions of population effects of hypoxia relies on the realism of the components of the model. While this approach can be very effective in isolating the effects of hypoxia, it is also time consuming as the data collection and modeling develops and evolves together.

Generality of Model Results
Our analysis is sound in some respects and, when trying to generalize the results, becomes more speculative in other aspects. We consider our estimation of the population effects on croaker in the NWGOM, when viewed as a relative change in the long-term averaged abundance, as being scientifically sound and defensible. We developed the components of the model in order to make such predictions and we erred on the side of modeling smaller, rather than larger, effects of low DO on individuals. We included density-dependent mortality to allow for long-term predictions, allowed for effective avoidance behavior, and developed growth, mortality, and fecundity effects submodels specifically for fluctuating exposures.
However, as we try to generalize our results, we must recognize caveats in the modeling. We used 3 years of DO output from FVCOM to represent intermediate and severe hypoxia years. More years that have different spatial and temporal patterns in hypoxia would further assure our results were robust. Similarly, the simulated DO in the croaker model still only varied over kilometers, while finer-scale variability could be important. We relied on a suite of models to generate DO and also changes in chlorophyll in response to reduced nutrient loadings; each carries their assumptions and limitations into and through our croaker modeling results. The link between chlorophyll and croaker food is also more complicated than our simple proportional relationships using multipliers. We are comfortable that we used available data and information to represent hypoxia and food effects sufficiently well to consider our predicted responses of croaker to be valid for historical conditions.
When we look more generally at the modeling, there are additional assumptions that cannot be addressed as easily. First, our analysis is a single-species approach that ignores indirect effects of hypoxia on prey and predators (i.e., the food web within which croaker live). This is a major assumption that is often used in fisheries management (Rose and Cowan 2003), but that has received criticism for oversimplifying the environmental conditions and food web interactions (e.g., Pikitch et al. 2004). Second, we used a species that is well studied and a major part of the fish community in Louisiana coastal waters (Monk et al. 2015), but also which will inhabit low DO waters that elicit growth and fecundity effects. The more we incorporated realism for croaker, the more difficult it becomes to extrapolate the results to other species (precision and accuracy versus generality, a la ' Levins 1966-see Odenbaugh 2003. How all these results can be extrapolated to other fish species with other avoidance triggers and behaviors (Gray et al. 2002;Wannamaker and Rice 2000) needs further evaluation and testing. In addition, while we allowed croaker to very effectively avoid low DO horizontally in our 2-D model, we did not allow for vertical movements. While there is evidence that croaker do not extensively move vertically as part of avoidance (Craig 2012), some other species are likely more inclined to use vertical movement to avoid low DO (Hazen et al. 2009). Simulating realistic 3-D behavioral movement of fish remains a challenge (LaBone 2016).
Our analysis showed the importance of realistically specifying the spatial and temporal variation in DO if our, or similar approaches, are further applied to additional species and systems. Higher variability in DO under FVCOM compared to idealized hypoxia emphasized the importance of specifying the neighborhood size underlying avoidance over the cue for selecting a destination cell and shifted the sensitivity of population responses from reproduction with little effects on mortality to mortality and reproduction. Careful specification of DO concentrations and how they relate to the resolution of the model grid, time step, and avoidance behavior and to the exposures used in laboratory experiments for growth, mortality, and reproduction effects is needed for each application. A more integrated set of models that couple hydrodynamics, water quality, zooplankton, and benthos would provide more consistent model inputs to fish population and food web models and enable more seamless comparisons of nutrient loadings and associated hypoxia effects on fish resources. Such quantification is needed for effective management (Bianchi et al. 2010b;Task Force 2013).