1 Introduction

Fish population (aka “stock”) assessment models determine the impact of past fishing on the historical and current abundance of the population, evaluate sustainable rates of removals (catch), and project future levels of catch reflecting one or more risk-averse catch rules. These catch rules are codified in regional Fishery Management Plans according to requirements of the Sustainable Fisheries Act. In the U.S., approximately 500 federally managed fish and shellfish populations are managed under approximately 50 Fishery Management Plans. About 200 of these populations are assessed each year, based on a prioritized schedule for their current status. Despite this, many minor species have never been quantitatively assessed. Although the pace is slower than that for weather forecasting, fish stock assessments are operational models for fisheries management.

Assessment models typically assimilate annual catches, data on fish abundance from diverse surveys and fishery sources, and biological information regarding fish body size and proportions at age. A suite of models is available depending on the degree of data availability and unique characteristics of the fish population or its fishery. Where feasible, environmental time series are used as indicators of changes in population or observation processes, especially to improve the accuracy of the projections of abundance and sustainable catch into the future. Such linkages are based principally on correlations given the challenge of conducting field observations on an appropriate scale. The frontier of model development is in the rapid estimation of parameters to include random temporal effects, in the simultaneous modeling of a suite of interacting species, and in the explicit treatment of the spatial distribution of the population.

Assessment models are loosely coupled to other models. For example, an ocean-temperature or circulation model or benthic-habitat map may be directly included in the pre-processing of the fish abundance survey. A time series of a derived ocean factor, like the North Atlantic Oscillation, can be included as an indicator of a change in a population process. Output of a multi-decadal time series of derived fish abundance can be an input to ecosystem and economic models to better understand cumulative impacts and benefits.

Stock Synthesis is an age- and size-structured assessment model in the class of models termed integrated analysis models. Stock Synthesis has evolved since its initial inception in order to model a wide range of fish populations and dynamics. The most recent major revision to Stock Synthesis occurred in 2016, when v.3.30 was introduced. This new version of Stock Synthesis required major revisions to the input files relative to earlier versions (see the Converting Files section for more information). The acronym for Stock Synthesis has evolved over time with earlier versions being referred to as SS2 (Stock Synthesis v.2.xx) and newer versions as SS3 (Stock Synthesis v.3.xx).

SS3 has a population sub-model that simulates a stock’s growth, maturity, fecundity, recruitment, movement, and mortality processes, an observation sub-model estimates expected values for various types of data, a statistical sub-model characterizes the data’s goodness of fit and obtains best-fitting parameters with associated variance, and a forecast sub-model projects needed management quantities. SS3 outputs the quantities, with confidence intervals, needed to implement risk-averse fishery control rules. The model is coded in C++ with parameter estimation enabled by automatic differentiation (admb). Windows, Linux, and iOS versions are available. Output processing and associated tools are in R, and a graphical interface is in QT. SS3 executables and support material is available on GitHub. The rich feature set in SS3 allows it to be configured for a wide range of situations. SS3 has become the basis for a large fraction of U.S. assessments and many other assessments around the world.

This manual provides a guide for using SS3. The guide contains a description of the input and output files and usage instructions. An overview and technical description of the model itself is in Methot and Wetzel (2013). However, SS3 has continued to evolve and grow since the publication in 2013, with this manual reflecting the most up-to-date information regarding SS3. The model and a graphical user interface are available on GitHub with older archived versions also available online at vlab. The vlab site also provides a user forum for posting questions and for accessing various additional materials. An output processor package, r4ss, in R is available for download from GitHub.

Additional guidance for new users can be found on the ss3-website which contains tutorials for getting started and building your own models as well as topic-focused vignettes.

To learn more about how to use Stock Synthesis, see the SS3 website for tutorials to get started and build your own models as well as topic-focused vignettes.

1.1 How To Cite

Please cite Stock Synthesis as:

Methot, R.D. and Wetzel, C.R. (2013). Stock Synthesis: A biological and statistical framework for fish stock assessment and fishery management. Fisheries Research, 142: 86-99. https://doi.org/10.1016/j.fishres.2012.10.012

Please cite the Stock Synthesis User Manual as:

Methot, R. D., Jr., C. R. Wetzel, I. G. Taylor, and K. Doering. (2020). Stock Synthesis User Manual Version 3.30.15. U.S. Department of Commerce, NOAA Processed Report NMFS-NWFSC-PR-2020-05. https://doi.org/10.25923/5wpn-qt71

2 File Organization

2.1 Input Files

1. starter.ss: required file containing filenames of the data file and the control file plus other run controls (required).

2. datafile: file containing model dimensions and the data (required)

3. control file: file containing set-up for the parameters (required)

4. forecast.ss: file containing specifications for reference points and forecasts (required)

5. ss3.par (previously ss.par): previously created parameter file that can be read to overwrite the initial parameter values in the control file (optional)

6. wtatage.ss: file containing empirical input of body weight by fleet and population and empirical fecundity-at-age (optional)

7. runnumber.ss: file containing a single number used as run number in output to CumReport.sso and in the processing of profilevalues.ss (optional)

8. profilevalues.ss: file contain special conditions for batch file processing (optional)

2.2 Output Files

1. data\_echo.ss\_new: Contains the input data as read by the model. In model versions prior to v.3.30.19 a single data.ss\_new file was created that included the echoed data, the expected data values (data\_expval.ss), and any bootstrap data files selected (data\_boot\_x.ss).

2. data\_expval.ss: Contains the expected data values given the model fit. This file is only created if the value for “Number of datafiles to produce” in the starter file is set to 2 or greater.

3. data\_boot\_x.ss: A new data file filled with bootstrap data based on the original input data and variances. This file is only created if the value in the “Number of datafiles to produce” in the starter file is set to 3 or greater. A separate bootstrap data file will be written for the number of bootstrap data file requests where x in the file name indicates the bootstrap simulation number (e.g., data\_boot\_001.ss, data\_boot\_002.ss, ...).

4. control.ss\_new: Updated version of the control file with final parameter values replacing the initial parameter values.

5. starter.ss\_new: New version of the starter file with annotations.

6. Forecast.ss\_new: New version of the forecast file with annotations.

7. warning.sso: This file contains a list of warnings generated during program execution. Starting in v.3.30.20 warnings are categorized into either “Note” or “Warning”. An item marked as “Note” denotes settings that the user may want to revise but do not require any additional changes for the model to run. Items marked with “Warning” are items that may or may not have allowed the model to finish running. Items with a fatal warning caused the model to fail during either reading input files or calculations. Warnings classified as error or adjustment may be causing calculation issues, even if the model was able to finish reading file and running, and should be addressed the user.

8. echoinput.sso: This file is produced while reading the input files and includes an annotated echo of the input. The sole purpose of this output file is debugging input errors.

9. Report.sso: This file is the primary report file.

10. ss\_summary.sso: Output file that contains all the likelihood components, parameters, derived quantities, total biomass, summary biomass, and catch. This file offers an abridged version of the report file that is useful for quick model evaluation. This file is only available in v.3.30.08.03 and greater.

11. CompReport.sso: Observed and expected composition data in a list-based format.

12. Forecast-report.sso: Output of management quantities and for forecasts.

13. CumReport.sso: This file contains a brief version of the run output, output is appended to current content of file so that the results of several runs can be collected together. This is useful when a batch of runs is being processed.

14. Covar.sso: This file replaces the standard admb ss.cor with an output of the parameter and derived quantity correlations in database format.

15. ss3.par (previously ss.par): This file contains all estimated and fixed parameters from the model run.

16. ss.std, ss.rep, ss.cor, etc.: Standard admb output files.

17. checkup.sso: Contains details of selectivity parameters and resulting vectors. This is written during the first call of the objective function.

18. Gradient.dat: New for v.3.30, this file shows parameter gradients at the end of the run.

19. rebuild.dat: Output formatted for direct input to Andre Punt’s rebuilding analysis package. Cumulative output is output to REBUILD.SS (useful when doing mcmc or profiles).

20. SIS\_table.sso: Output formatted for reading into the nmfs sis.

21. Parmtrace.sso: Parameter values at each iteration.

22. posteriors.sso, derived\_posteriors.sso, posterior\_vectors.sso: Files associated with mcmc.

3 Starting Stock Synthesis

SS3 is typically run through the command line interface, although it can also be called from another program, R, the ssi, or a script file (such as a DOS batch file). SS3 is compiled for Windows, Mac, and Linux operating systems. The memory requirements depend on the complexity of the model you run, but in general, SS3 will run much slower on computers with inadequate memory. See Running Stock Synthesis for additional notes on methods of running SS3.

Communication with the program is through text files. When the program first starts, it reads the file starter.ss, which typically must be located in the same directory from which SS3 is being run. The file starter.ss contains required input information plus references to other required input files, as described in the File Organization section. The names of the control and data files must match the names specified in the starter.ss file. File names, including starter.ss, are case-sensitive on Linux and Mac systems but not on Windows. The echoinput.sso file outputs how the executable reads each input file and can be used for troubleshooting when trying to set up a model correctly. Output from SS3 consists of text files containing specific keywords. Output processing programs, such as Excel or R, can search for these keywords and parse the specific information located below that keyword in the text file.

4 Converting Files from Stock Synthesis v.3.24

Converting files from version 3.24 to version 3.30 can be performed by using the program ss_trans.exe. This executable takes v.3.24 files as input and will output v.3.30 input and output files. SS_trans executables are available for v.3.30.01 - v.3.30.17. The transitional executable was phased out with v.3.30.18. If a model needs to be converted from v.3.24 to a recent version, one should use the v.3.30.17 ss_trans.exe available from the v.3.30.17 release page on GitHub to convert the files and then any additional adjustments needed between v.3.30.17 and newer versions should be done by hand. To see the changes that need to be made between v.3.30.17 and the latest release of SS3, please see the change log for v.3.30.19 onward as well as the Excel version of the change log for versions prior to v.3.30.19.

The following file structure and steps are recommended for converting model files:

1. Create “transition” folder. Place the 4 main model files (control, data, starter, and forecast) from version SS3 v.3.24 within the transition folder along with the SS3 transition executable (ss_trans.exe). One tip is to use the control.ss_new from the SS3 v.3.24 estimated model rather than the original control file which will set all parameter values at the previous estimated mle parameters. Run the transition executable with a command like ss_trans -stopph -1 which will write ss_new files but not go through the population dynamics or produce other output (whereas -stopph 0 will run through the dynamics once without estimation and produce Report.sso and other output files, but may not produce the ss_new files if there is any issue with the model setup).

2. Create “converted” folder. Place the ss_new (data.ss_new, control.ss_new, starter.ss_new, forecast.ss_new) files created by the transition executable contained within the “transition” folder into this new folder. Rename the ss_new files to the appropriate suffixes and change the names in the starter.ss file accordingly.

3. Review the control (control.ss_new) file to determine that all model functions converted correctly. The structural changes and assumptions for some advanced model features are too complicated to convert automatically. See below for some known features that may not convert. When needed, it is recommended to modify the control.ss_new file, the converted control file, for only the features that failed to convert properly.

4. Change the max phase to a value greater than the last phase in which the parameter is set to be estimated within the control file. Run the new v.3.30 executable (ss3.exe) within the “converted” folder using the renamed ss_new files created from the transition executable.

5. Compare likelihood and model estimates between the v.3.24 and v.3.30 model versions.

6. If desired, update to versions of Stock Synthesis > v.3.30.17 by running the new v.3.30 input files with the higher executable.

Some options have been substantially changed in v.3.30, which impedes the automatic converting of v.3.24 model files. Known examples of v.3.24 options that cannot be converted, but for which better alternatives are available in v.3.30 are:

1. The use of catchability (Q) deviations,

2. Complex birth seasons,

3. Environmental effects on spawner-recruitment parameters,

4. Setup of time-varying quantities for models that used the no-longer-available features (e.g., logistic bound constraint).

5 Starter File

5.1 Reading the Manual’s format

SS3 begins by reading the file starter.ss. The starter file contains need information on the names of the control and data files, run conditions, and output specifications. The term COND appears in the “Typical Value” column of this documentation (it does not actually appear in the model files), it indicates that the following section is omitted except under certain conditions, or that the factors included in the following section depend upon certain conditions. In most cases, the description in the definition column is the same as the label output to the ss_new files.

5.2 Terminology for Fishing Mortality, \(F\)

Here we introduce some terminology related to fishing mortality, \(F\). This will provide context for some of the quantities that will be read from the starter file and used throughout the document

\(f\) is fleet.

\(t\) is a time step; continuous across years \(y\) and seasons \(s\); equivalent to year if only 1 season.

\(a\) is age.

\(s_{t,f,a}\) is age-specific selectivity for a fleet. If selectivity is length-specific, then age-specific selectivity due to length-selectivity is calculated as the dot product across length bins of length selectivity and the normal (or log-normal) distribution of length-at-age. If selectivity is both length- and age-based, which is an entirely normal concept in SS3, then age selectivity due to length selectivity is calculated first, then multiplied by the direct age selectivity. This compound age selectivity is used in the mortality calculations and is reported as Asel2 in report:32 of Report.sso. Selectivity can be sex-specific, and different growth morphs and platoons can have different age-selectivity due to the effect of length-selectivity on their unique size-at-age. This added dimension, \(g\), for biological group is not included in the nomenclature here but exists in all the SS3 calculations.

\(F_{t,f,a}\) is fishing mortality at age for fleet \(f\). There is no subscript for area because each fleet is defined to operate in only one area.

\(F_{t,f}'\) is a fleet’s fishing mortality for the age that has selectivity equal to 1.0. This is also termed F’ or \(\text{full\_}F\) in the SS3 system. If your model is using parameters for \(F\), then the parameter values are for the \(F'\). Note that some selectivity curves, like double normal, are explicit about having a maximum of 1.0. But other curves like logistic and combinations of length-selectivity and growth, may produce an age-selectivity curve that never reaches 1.0 and time-varying non-parametric selectivity will produce values > 1.0 routinely. In all cases, the resultant \(F_{t,f,a}\) comes from \(F_{t,f}' * s_{t,f,a}\), so the range of the \(F'\) compensates for the scale of the \(s\).

Apical selectivity is the maximum age-specific selectivity and is not explicit in any internal calculation in SS3, it is just for reporting. If selectivity has a maximum value of 1.0, then \(\text{apical\_}F\) and \(\text{full\_}F\) are identical.

Fully-selected age range is not explicitly used in SS3, especially because SS3 applications routinely have multiple fleets with different selectivity patterns that may have little overlap.

Fbar is the average \(F\) over a user specified range of ages, implicitly the fully-selected range for the total \(F\) from all the fleets. Some SS3 output options will display Fbar.

\(\text{Annual\_}F\) is essentially the same as Fbar and is an output quantity.

\(F\text{\_std}\) is an output quantity that may be based on \(\text{annual\_}F\) or other calculated quantities like exploitation rate. Importantly, the output values of \(F\text{\_std}\) may be presented as a ratio relative to an equivalent benchmark (reference point) quantity; e.g., \(F / F_{MSY}\). Further, the variance of \(F\text{\_std}\) will be calculated and output.

\(C_{t,f}\) is fleet-specific catch in a time step.

\(B_{t,f}\) is fleet specific available biomass, e.g., total biomass filtered by fleet-specific age selectivity, \(s_{t,f,a}\). Note that this is not adjusted by the \(max(s_{t,f,a})\).

6 Forecast File

The specification of options for forecasts is contained in the mandatory input file named forecast.ss. See Forecast Module: Benchmark and Forecasting Calculations for additional details.

The term COND appears in the “Typical Value” column of this documentation (it does not actually appear in the model files) and indicates that the following section is omitted except under certain conditions, or that the factors included in the following section depend upon certain conditions. In most cases, the description in the definition column is the same as the label output to the ss_new files.

6.2 Including a New Fleet in the Forecast

As of v.3.30.16 users can have a forecast fleet without catches during the modeled period. Previously, fleets in the forecast period were required to have input catches at some amount during the modeled period. SS3 now has capability to have a fleet with no input catches during the modeled period that could be used as a fleet during the forecast.

6.3 Benchmark Calculations

This feature of SS3 is designed to calculate an equilibrium fishing rate intended to serve as a proxy for the fishing rate that would provide maximum sustainable yield (\(F_{MSY}\)). Then in the forecast module these fishing rates can be used in the projections.

Four reference points can be calculated by SS3. The first is the estimate of \(F_{MSY}\) within the model, while the other reference points use proxies or an alternative estimated point.

• \(F_{MSY}\): Search for the \(F\) that produces maximum equilibrium (e.g., dead catch).

• \(F_{SPR}\): Search for the \(F\) that produces spawning biomass per recruit this is a specific fraction, termed \(SPR_{TARGET}\), of spawning biomass per recruit under unfished conditions. Note that this is in relative terms so it does not take into account the spawner-recruit relationship.

• \(F_{B_{TARGET}}\): Search for the \(F\) that produces an absolute spawning biomass that is a specified fraction, termed relative biomass target, of the unfished spawning biomass. Note that this is in absolute terms so takes into account the spawner-recruit relationship.

• \(F_{0.10}\): Search for the \(F\) that produces a slope in yield per recruit, \(dY/dF\), that is 10% of the slope at the origin. Note that with SS3, this option is mutually exclusive with \(F_{B_{TARGET}}\). Only one will be calculated and the one that is calculated can serve as the proxy for \(F_{MSY}\) and forecasting. The \(F_{0.10}\) search can fail if bycatch fleets are used and the bycatch setting includes the fleet’s catch in the catch to be optimized.

6.3.0.1 Estimation

Each of the potential reference points is calculated by searching across a range of \(F\) multiplier levels, calculating equilibrium biomass and catch at that \(F\), using Newton-Raphson method to calculate a better \(F\) multiplier value, and iterating a fixed number of times to achieve convergence on the desired level.

6.3.0.2 Calculations

The calculation of equilibrium biomass and catch uses the same code that is used to calculate the virgin conditions and the initial equilibrium conditions. This equilibrium calculation code takes into account all morph, timing, biology, selectivity, and movement conditions as they apply while doing the time series calculations. You can verify this by running SS3 to calculate \(F_{MSY}\) then hard-wire initial \(F\) to equal this value, use the F_method approach 2 so each annual \(F\) is equal to \(F_{MSY}\) and then set forecast \(F\) to be the same \(F_{MSY}\). Then run SS3 without estimation and no recruitment deviations. You should see that the population has an initial equilibrium abundance equal to \(B_{MSY}\) and stays at this level during the time series and forecast.

6.3.0.3 Catch Units

For each fleet, SS3 always calculates catch in terms of biomass (mt) and numbers (1000s) for encountered (selected) catch, dead catch, and retained catch. These three categories differ only when some fleets have discarding or are designated as a bycatch fleet. SS3 uses total dead catch biomass as the quantity that is principally reported and the quantity that is optimized when searching for \(F_{MSY}\). The quantity “dead catch” may occasionally be referred to as “yield”.

6.3.0.4 Biomass Units

The principle measure of fish abundance, for the purpose of reference point calculation, is female reproductive output. This is referred to as ssb and sometimes just “B” because the typical user settings have one unit of reproductive output (fecundity) per kg of mature female biomass. So when the output label says \(B_{MSY}\), this is actually the female reproductive output at the proxy for \(F_{MSY}\).

6.3.0.5 Fleet Allocation

An important concept for the reference point calculation is the allocation of fishing rate among fleets. Internally, this is benchmark years relative \(F\) (\(text{Bmark\_rel}F\) (\(f,s\))) and it is the fraction of the \(F\) multiplier assigned to each fleet, \(f\) and season, \(s\). The value, \(F\text{mult} * \text{Bmark\_rel}F\)(\(f,s\)), is the \(F\) level for a particular fleet in a particular season and for the age that has a selectivity of 1.0. Other ages will have different \(F\) values according to their selectivity.

• The \(\text{Bmark\_rel}F\) values can be calculated by SS3 from a range of years specified in the input for Benchmark Years or it can be set to be the same as the Forecast_RelF, which in turn can be based on a range of years or can be input as a set of fixed values.

• The biology years selected for the \(\text{Bmark\_rel}F\) calculations can have an effect on the standard deviation estimated, even in models with no time-varying biology. It is recommended to use the first model year except when biology is time-varying. Note that when biology is time-varying, the fecundity vector, which is used to calculate ssb, will be updated in every iteration for every year that has time-varying biology.

• Note for Bycatch Fleets: The forecast \(F\) for a bycatch fleet is normally set in the bycatch setup in the data file to be the average of a range of years or a fixed input value, but year-specific \(F\) values can be set at the end of the forecast file in the Forecast catch input section. If this is done, that \(F\) value is not adjusted by changes to the \(F\) multiplier. This allows the user to treat a bycatch fleet as a constant background \(F\) while the optimal \(F\) for other fleets is sought. Also, for bycatch fleets, there is user control for whether the dead catch from the bycatch fleet is included in the total dead catch that is optimized when searching for \(F_{MSY}\).

6.3.0.6 Virgin vs. Unfished Spawning Biomass

The concept of unfished spawning biomass, (written as SSB_unfished in SS3 input and output files), is important to the reference points calculations. Unfished spawning biomass can be potentially different from virgin spawning biomass (written as SSB_virgin in SS3 output files).

• Virgin spawning biomass is calculated from the parameter values associated with the start year of the model configuration and it serves as the basis from which the population model starts and the basis for calculation of stock depletion.

• Unfished spawning biomass can be calculated for any year or range of years, so can change over time as \(R_{0}\), steepness, or biological parameters change.

• In the reference points calculation, the Benchmark Years input specifies the range of time over which the mean of various quantities are taken to calculate the reference points. For biology, selectivity, \(F\)s, and movement the mean values are the year-specific derived quantities. But for the stock-recruitment parameters (\(R_{0}\) and steepness), the mean of the parameter values themselves is calculated over time.

• During the time series or forecast, the current year’s unfished spawning output is used as the basis for the spawner-recruitment curve against which deviations from the spawner-recruitment curve are applied. So, if \(R_{0}\) is made time-varying, then the spawner-recruit curve itself is changed. However, if the regime shift parameter is time-varying, then this is an offset from the spawner-recruitment curve and not a change in the curve itself. Changes in \(R_{0}\) will change year-specific reference points and change the expected value for annual recruitments, but changes in regime shift parameter only change the expected value for annual recruitments.

• In reporting the time series of depletion level, the denominator can be based on virgin spawning output or \(B_{MSY}\). Note that \(B_{MSY}\) is based on unfished spawning output for the specified range of Benchmark years, not on virgin spawning biomass.

6.4 Forecast Recruitment Adjustment

Recruitment during the forecast years sometimes needs to be set at a level other than that determined by the spawner-recruitment curve. One way to do this is by an environmental or block effect on the regime shift parameter. A more straightforward approach is now provided by the special forecast recruitment feature described here. There are 4 options provided for this feature. These are:

• 0 = Do nothing: this is the default and will invoke no special treatment for the forecast recruitments.

• 1 = Multiplier on spawner-recruitment: the expected recruitment from the srr is multiplied by this factor.

  • This is a multiplier, so null effect comes from a value of 1.0.

  • The order of operations is to apply the srr, then the regime effect, then this special forecast effect, then bias adjustment, then the deviations.

  • In the spawner recruit output of the Report.sso there are 4 recruitment values stored.

• 2 = Multiplier on virgin recruitment: The virgin recruitment is multiplied by this factor.

  • This is a multiplier, so null effect comes from a value of 1.0.

  • The order of operations is to apply any environmental or block effects to \(R_{0}\), then apply the special forecast effect, then bias adjustment, then the deviations.

  • Note that environmental or block effects on \(R_{0}\) are rare and are different from environment or block effects on the regime parameter.

• 3 = Mean recent recruitment: calculate the mean recruitment and use it during the forecast period.

  • Note that bias adjustment is not applied to this mean because the values going into the mean have already been bias adjusted.

This feature affects the expected recruitment in all years after the last year of the main recruitment deviations. This means that if the last year of main recruitment deviations is before end year, then the last few recruitments, termed “late”, are also affected by this forecast option. For example, option 3 would allow you to set the last 2 years of the time series and all forecast years to have recruitment equal to the mean recruitment for the last 10 years of the main recruitment era.

7 Data File

7.1 Overview of Data File

1. Dimensions (years, ages, number of fleets, number of surveys, etc.)

2. Fleet and survey names, timing, etc.

3. Catch data (biomass or numbers)

4. Discard totals or rate data

5. Mean body weight or mean body length data

6. Length composition set-up

7. Length composition data

8. Age composition set-up

9. Age imprecision definitions

10. Age composition data

11. Mean length-at-age or mean bodyweight-at-age data

12. Generalized size composition (e.g., weight frequency) data

13. Environmental data

14. Tag-recapture data

15. Stock composition (e.g., morphs identified by otolith microchemistry) data

16. Selectivity observations (new placeholder, not yet implemented)

7.2 Units of Measure

The normal units of measure are as follows:

• Catch biomass - metric tons

• Body weight - kg

• Body length - usually in cm, weight-at-length parameters must correspond to the units of body length and body weight

• Survey abundance - any units if catchability (Q) is freely scaled; metric tons or thousands of fish if Q has a quantitative interpretation

• Output biomass - metric tons

• Numbers - thousands of fish, because catch is in metric tons and body weight is in kg

• Spawning biomass (ssb)- metric tons of mature females if eggs/kg = 1 for all weights; otherwise has units that are based on the user-specified fecundity

7.3 Time Units

• Spawning:

  • Happens once per year at a specified date (in real months, 1.0 - 12.99). To create multiple spawning events per year, change the definition of a year (e.g., call a 6-month period a “year” when there are two spawning events about 6 months apart). However, revising the definition of year affects assignment of fish age, so it should not be used if age data are included.

• Recruitment:

  • Occurs at specified recruitment events that occur at user-specified dates (in real months, 1.0 - 12.99).

  • There can be one to many recruitment events across a year; each producing a platoon as a portion of the total recruitment.

  • A settlement platoon enters the model at age 0 if settlement is between the time of spawning and the end of the year; it enters at age 1 if settlement is after the first of the year; these ages at settlement can be overridden in the settlement setup.

• Timing

  • All fish advance to the next older integer age on January 1, no matter when they were born during the year. Consult with your ageing lab to assure consistent interpretation.

• Parameters:

  • Time-varying parameters are allowed to change annually, not seasonally.

  • Rates like growth and mortality are per year.

7.3.1 Seasons

Seasonal quantities in the model are calculated and treated in the following methods:

• Seasons are the time step during which constant rates apply.

• Catch and discard amounts are per season and \(F\) is calculated per season.

• The year can have just one annual season, or be subdivided into seasons of unequal length.

• Season duration is input in real months (1.0 - 12.99) and is converted into fractions of an annum. Annual rate values are multiplied by the per annum season duration.

• If the sum of the input season durations do not equal 12.0, then the input durations are divided by the total duration input to rescale season duration to equal 1.0.

• Allows for a special situation in which the year could be only 0.25 in duration (e.g., seasons as years) so that spawning and time-varying parameters can occur more frequently. Other durations are also possible (e.g., where a model year is really only 6 months). Note that real month inputs should always have the range 1.0 - 12.99, even in cases where a model year is not a true year. See the continuous seasonal recruitment section for more information on how to set up models where the model year duration is different from a true year.

7.3.2 Subseasons and Timing of Events

The treatment of subseasons in SS3 provide more precision in the timing of events compared to earlier model versions. In early versions, v.3.24 and before, there was effectively only two subseasons per season because the alk for each observation used the mid-season mean length-at-age and spawning occurred at the beginning of a specified season.

Time steps can be broken into subseason and the alk can be calculated multiple times over the course of a year:

• Even number (min = 2) of subseasons per season (regardless of season duration):

  • Two subseasons will mimic v.3.24

  • Specifying more subseasons will give finer temporal resolution, but will slow the model down, the effect of which is mitigated by only calculating growth as needed.

• Survey timing is now cruise-specific and specified in units of months (e.g., April 15 = 4.5; possible inputs are 1.0 to 12.99).

  • ss_trans.exe will convert year, season in v.3.24 format to year, real month in v.3.30 format.

• Survey integer season and spawn integer season assigned at run time based on real month and season duration(s).

• Growth and the alk is calculated at the beginning and mid-season or when there is a defined subseason with data observation.

• Fishery body weight uses mid-subseason growth.

• Survey body weight and size composition is calculated using the nearest subseason.

• Reproductive output now has specified spawn timing (in months fraction) and interpolates growth to that timing.

• Survey numbers calculated at cruise survey timing using \(e^{-z}\).

• Continuous \(Z\) for entire season. Same as applied in version v.3.24.

7.4 Terminology

The term COND appears in the “Typical Value” column of this documentation (it does not actually appear in the model files), it indicates that the following section is omitted except under certain conditions, or that the factors included in the following section depend upon certain conditions. In most cases, the description in the definition column is the same as the label output to the ss_new files.

7.5 Model Dimensions

7.6 Fleet Definitions

The catch data input has been modified to improve the user flexibility to add/subtract fishing and survey fleets to a model set-up. The fleet setup input is transposed so each fleet is now a row. Previous versions (v.3.24 and earlier) required that fishing fleets be listed first followed by survey only fleets. In SS3 all fleets have the same status within the model structure and each has a specified fleet type (except for models that use tag recapture data, this will be corrected in future versions). Available types are; catch fleet, bycatch only fleet, or survey.

7.6.0.1 Fleet Type

Define the fleet type (e.g., fishery fleet, survey fleet):

• 1 = fleet with input catches;

• 2 = bycatch fleet (all catch discarded) and invoke extra input for treatment in equilibrium and forecast;

• 3 = survey: assumes no catch removals even if associated catches are specified below. If you would like to remove survey catch set fleet type to option = 1 with specific month timing for removals (defined below in the “Timing” section); and

• 4 = predator (M2) fleet that adds additional mortality without a fleet \(F\) (added in v.3.30.18). Ideal for modeling large mortality events such as fish kills or red tide. Requires additional long parameter lines for a second mortality component (M2) in the control file after the natural mortality/growth parameter lines (entered immediately after the fraction female parameter line).

7.6.0.2 Timing

Timing for data observations:

• Fishery options:

  • -1: catch is treated as if it occurred over the whole season. SS3 may change the cpue data to occur in the middle of the season if it is specified otherwise (i.e., the cpue observations may have a different month in the data_echo.ss_new file). A user can override this assumption for specific data observations (e.g., length or age) by specifying a month. This option works well for fisheries where fishing is spread throughout the year.

  • 1: The fleet timing is not used and only the month value associated with each observation is relevant. This option works well for pulse fisheries that occurs over a small subset of months.

• Survey option, 1: The fleet timing is not used and only the month value associated with each observation is relevant (e.g., month specification in the indices of abundance or the month for composition data). This input should always be used for surveys.

7.6.0.3 Area

An integer value indicating the area in which a fleet operates.

7.6.0.4 Catch Units

Ignored for survey fleets, their units are read later:

• 1 = biomass (in metric tons); and

• 2 = numbers (thousands of fish).

See Units of Measure for more information.

7.6.0.5 Catch Multiplier

Invokes use of a catch multiplier, which is then entered as a parameter in the mortality-growth parameter section. The estimated value or fixed value of the catch multiplier is used to adjust the observed catch:

• 0 = No catch multiplier used; and

• 1 = Apply a catch multiplier which is defined as an estimable parameter in the control file after the cohort growth deviation in the biology parameter section. The model’s estimated retained catch will be multiplied by this factor before being compared to the observed retained catch.

A catch multiplier can be useful when trying to explore historical unrecorded catches or ongoing illegal and unregulated catches. The catch multiplier is a full parameter line in the control file and has the ability to be time-varying.

7.7 Bycatch Fleets

The option to include bycatch fleets was introduced in v.3.30.10. This is an optional input and if no bycatch is to be included in to the catches this section can be ignored.

A fishing fleet is designated as a bycatch fleet by indicating that its fleet type is 2. A bycatch fleet creates a fishing mortality, same as a fleet of type 1, but a bycatch fleet has all catch discarded, so the input value for retained catch is ignored. However, an input value for retained catch is still needed to indicate that the bycatch fleet was active in that year and season. A catch multiplier cannot be used with bycatch fleets because catch multiplier works on retained catch. SS3 will expect that the retention function for this fleet will be set in the selectivity section to type 3, indicating that all selected catch is discarded dead. It is necessary to specify a selectivity pattern for the bycatch fleet and, due to generally lack of data, to externally derive values for the parameters of this selectivity.

All catch from a bycatch fleet is discarded, so one option to use a discard fleet is to enter annual values for the amount (not proportion) that is discarded in each time step. However, it is uncommon to have such data for all years. An alternative approach that has been used principally in the U.S. Gulf of Mexico is to input a time series of effort data for this fleet in the survey section (e.g., effort is a “survey” of \(F\), for example, the shrimp trawl fleet in the Gulf of Mexico catches and discards small finfish and an effort time series is available for this fleet) and to input in the discard data section an observation for the average discard over time using the super year approach. Another use of bycatch fleet is to use it to estimate effect of an external source of mortality, such as a red tide event. In this usage there may be no data on the magnitude of the discards and SS3 will then rely solely on the contrast in other data to attempt to estimate the magnitude of the red tide kill that occurred. The benefit of doing this as a bycatch fleet, and not a block on natural mortality, is that the selectivity of the effect can be specified.

Bycatch fleets are not expected to be under the same type of fishery management controls as the retained catch fleets included in the model. This means that when SS3 enters into the reference point equilibrium calculations, it would be incorrect to have SS3 re-scale the magnitude of the \(F\) for the bycatch fleet as it searches for the \(F\) that produces, for example, F35%. Related issues apply to the forecast. Consequently, a separate set of controls is provided for bycatch fleets (defined below). Input is required for each fleet designated as fleet type = 2.

If a fleet above was set as a bycatch fleet (fleet type = 2), the following line is required:

The above example set-up defines one fleet (fleet number 2) as a bycatch fleet with the dead catch from this fleet to not be included in the search for msy (b: Include in msy = 2). The level of \(F\) from the bycatch fleet in reference point and forecast is set to the mean (c: \(F\text{mult}\) = 3) of the estimated \(F\) for the range of years from 1982-2010.

7.7.0.1 Fleet Index

Fleet number for which to include bycatch catch. Fleet number is assigned within the model based on the order of listed fleets in the Fleet Definition section. If there are multiple bycatch fleets, then a line for each fleet is required in the bycatch section.

7.7.0.2 Include in msy

The options are:

• 1 = dead fish in msy, abc, and other benchmark and forecast output; and

• 2 = omit from msy and abc (but still include the mortality).

7.7.0.3 \(F\) Multiplier (\(F\text{mult}\))

The options are:

• 1 = \(F\) multiplier scales with other fleets;

• 2 = bycatch \(F\) constant at input value in column d; and

• 3 = bycatch \(F\) from range of years input in columns d and e.

7.7.0.4 \(F\) or First Year

The specified \(F\) or first year for the bycatch fleet.

7.7.0.5 \(F\) or Last Year

The specified \(F\) or last year for the bycatch fleet.

7.7.0.6 Not Used

This column is not yet used and is reserved for future features.

7.7.0.7 Bycatch Fleet Usage Instructions and Warnings

When implementing a bycatch fleet, changes to both the data and control file are needed.

The needed changes to the data file are:

1. Fleet type - set to value of 2.

2. Set bycatch fleet controls per information above.

3. Catch input - you must enter a positive value for catch in each year/season that you want a bycatch calculated. The entered value of catch will be ignored by SS3, it is just a placeholder to invoke creating an \(F\).

   1. Initial equilibrium - you may want to enter the bycatch amount as retained catch for the initial equilibrium year because there is no option to enter initial equilibrium discard in the discard section.

4. Discard input - It is recommended to enter the amount of discard to assist SS3 in estimating the \(F\) for the bycatch fleet.

5. Survey input - It is useful, but not absolutely necessary, to enter the effort time series by the bycatch fleet to assist SS3 in estimating the annual changes in \(F\) for the bycatch fleet.

The needed changes to the control file are:

1. The \(F\) method must be set to 2 in order for SS3 to estimate \(F\) with having information on retained catch.

2. Selectivity -

   1. A selectivity pattern must be specified and fixed (or estimated if composition data is provided).

   2. The discard column of selectivity input must be set to a value of 3 to cause all catch to be discarded.

In v.3.30.14 it was identified that there can be an interaction between the use of bycatch fleets and the search for the \(F_{0.1}\) reference point which may results in the search failing. Changes to the search feature were implemented to make the search more robust, however, issue may still be encountered. In these instances it is recommended to not select the \(F_{0.1}\) reference point calculation in the forecast file.

7.8 Predator Fleets

Introduced in v.3.30.18, a predator fleet provides the capability to define an entity as a predator that adds additional mortality (\(M2\), i.e., the predation mortality) to the base natural mortality. This new capability means that previous use of bycatch fleets to mimic predators (or fish kills, e.g., due to red tide) will no longer be necessary. The problem with using a bycatch fleet as a predator was that it still created an \(F\) that was included in the reporting of total \(F\) even if the bycatch was not included in the msy search.

For each fleet that is designated as a predator, a new parameter line is created in the mg parameter section in the control file. This parameter will have the label M2_pred1, where the “1” is the index for the predator (not the index of the fleet being used as a predator). More than one predator can be included. If the model has > 1 season, it is normal to expect \(M2\) to vary seasonally. Therefore, only if the number of seasons is greater than 1, follow each \(M2\) parameter with number of season parameters to provide the seasonal multipliers. These are simple multipliers times \(M2\), so at least one of these needs to have a non-estimated value. The set of multipliers can be used to set \(M2\) to only operate in one season if desired. If there is more than one predator fleet, each will have its own seasonal multipliers. If there is only 1 season in the model, then no multiplier lines are included.

Three types of data relevant to \(M2\) can be input:

• Total kill (as discard in the data file): \(M2\) is a component of \(Z\), so \(M2/Z\) can be used to calculate the amount of the total kill that is attributable to \(M2\). This is completely analogous to calculating catch for the fishing fleets. The total kill (e.g., consumption) is output to the discard array. If data on the total kill by the \(M2\) predator is available, it can be input as observed “discard” for this fleet and thus included in the total log likelihood to estimate the magnitude of the \(M2\) parameter.

• Predator effort (as a survey index in the data file): \(M2\) is a rate analogous to \(F\), so the survey of \(F\) approach (survey units = 2) can be used to input predator abundance as an indicator of the “effort” that produced the \(M2\). Like all surveys, this survey of \(M2\) will also need a Q specification. Note that in the future we can explore improved options for this Q.

• Predated age-length composition (as length or age composition data in the data file): \(M2\) “eats” the modeled fish, so gut contents or other sources may have size and/or age composition data which may be input to estimate selectivity of the \(M2\) source.

With the input of data on the time series of total kill or predator effort, it should be possible to estimate annual deviations around the base \(M2\) for years with data. If the \(M2\) time series is instead driven by environmental data, then also including data on kill or effort can provide a means to view consistency between the environmental time series and the additional data sets. Output of \(M2\) is found in a Report.sso section labeled predator (\(M2\)). In the example below, the \(M2\) seasonal multiplier was defined to have random deviations by year. This allowed multipliers plus \(M2\) itself to closely match the input consumption amounts (288 mt of consumption per season, the fit can be examined by looking at the discard output report).

7.9 Catch

After reading the fleet-specific indicators, a list of catch values by fleet and season are read in by the model. The format for the catches is year and season that the catch is attributed to, fleet, a catch value, and a year-specific catch standard error. Only positive catches need to be entered, so there is no need for records corresponding to all years and fleets. To include an equilibrium catch value for a fleet and season, the year should be noted as -999. For each non-zero equilibrium catch value included, a short parameter line is required in the initial \(F\) section of the control file.

There is no longer a need to specify the number of records to be read; instead the list is terminated by entering a record with the value of -9999 in the year field. The updated list based approach extends throughout the data file (e.g., catch, length- and age-composition data), the control file (e.g., lambdas), and the forecast file (e.g., total catch by fleet, total catch by area, allocation groups, forecasted catch).

In addition, it is possible to collapse the number of seasons. So, if a season value is greater than the number of seasons for a particular model, that catch is added to the catch for the final season. This is one way to easily collapse a seasonal model into an annual model. The alternative option is to the use of season = 0. This will cause SS3 to distribute the input value of catch equally among the number of seasons. SS3 assumes that catch occurs continuously over seasons and hence is not specified as month in the catch data section. However, all other data types will need to be specified by month.

The format for a 2 season model with 2 fisheries looks like the table below. Example is sorted by fleet, but the sort order does not matter. In data.ss_new, the sort order is fleet, year, season.

• Catch can be in terms of biomass or numbers for each fleet, but cannot be mixed within a fleet.

• Catch is retained catch (aka landings). If there is discard also, then it is handled in the discard section below. This is the recommended setup which results in a model estimated retention curve based upon the discard data (specifically discard composition data). However, there may be instances where the data do not support estimation of retention curves. In these instances catches can be specified as all dead (retained + discard estimates).

• If there are challenges to estimating discards within the model, catches can be input as total dead without the use of discard data and retention curves.

• If there is reason to believe that the retained catch values underestimate the true catch, then it is possible in the retention parameter set up to create the ability for the model to estimate the degree of unrecorded catch. However, this is better handled with the new catch multiplier option.

7.10 Surveys and Indices

Indices are data that are compared to aggregate quantities in the model. Typically, the index is a measure of selected fish abundance, but this data section also allows for the index to be related to a fishing fleet’s \(F\), or to another quantity estimated by the model. The first section of the “Indices” setup contains the fleet number, units, error distribution, and whether additional output (sd Report) will be written to the Report file for each fleet that has index data.

7.10.0.1 Units

The options for units for input data are:

• 0 = numbers;

• 1 = biomass;

• 2 = \(F\); and

  • Note the \(F\) option can only be used for a fishing fleet and not for a survey, even if the survey selectivity is mirrored to a fishing fleet. The values of these effort data are interpreted as proportional to the level of the fishery \(F\) values. No adjustment is made for differentiating between continuous \(F\) values versus exploitation rate values coming from Pope’s approximation. A normal error structure is recommended so that the input effort data are compared directly to the model’s calculated \(F\), rather than to \(ln(F)\). The resultant proportionality constant has units of 1/Q where Q is the catchability coefficient. For more information see the section on Predator effort.

• \(>=\) 30 = Special survey types. These options bypass the calculation of survey selectivity so the no selectivity parameter are required and age/length selectivity pattern should be set as 0. A catchability parameter line in the control file will be required for each special survey. Special survey types 31, 32, and 36 relate to recruitment deviations. Before v.3.30.22, the expected value for observations before recdev_start or after recdev_end were null. With v.3.30.22, expected values are now based on recruitment deviations for all years and suggestions are included in warnings.sso if observations occur outside the range of active recruitment deviations. The expected values for these types are:

  • 30 = spawning biomass/output (e.g., for an egg and larvae survey);

  • 31 = exp(recruitment deviation), useful for environmental index affecting recruitment;

  • 32 = spawning biomass * exp(recruitment deviation), for a pre-recruit survey occurring before density-dependence;

  • 33 = recruitment, age-0 recruits;

  • 34 = depletion (spawning biomass/virgin spawning biomass);

    • Special survey option 34 automatically adjusts phases of parameters. To use the depletion survey approach, the user will need to make the following revisions to the SS3 data file: 1) add a new survey fleet, 2) define the survey type as option 34, 3) add two depletion survey data points, and initial unfished set equal to 1 for an unfished modeled year and one for a later year with the depletion estimates, 4) set the input cv value for each survey data point to a low value (e.g., 0.0001) to force the model to fit these data, and in the control file 5) add the survey to the control file in the Q set-up and selectivity sections with float set to 0 with parameter value set to 0.

    • There are options for additional control over this in the control file catchability setup section under the link information bullet where:

      • 0 = add 1 to phases of all parameters. Only \(R_{0}\) active in new phase 1. Mimics the default option of previous model versions;

      • 1 = only \(R_{0}\) active in phase 1. Then finish with no other parameters becoming active; useful for data-limited draws of other fixed parameters. Essentially, this option allows SS3 to mimic dbsra; and

      • 2 = no phase adjustments, can be used when profiling on fixed \(R_{0}\).

    • Warning: the depletion survey approach has not been tested on multiple area models. This approach may present challenges depending upon the dynamics within each area.

  • 35 = survey of a deviation vector (\(e(survey(y)) = f(parm\_dev(k,y))\)), can be used for an environmental time series that serves as an index for a parameter deviation vector. The selected deviation vector is specified in Q section of the control file. The index of the deviation vector to which the index is related is specified in the 2nd column of the Q setup table (see Catchability);

  • 36 = recruitment deviation

7.10.0.2 Error Distribution

The options for error distribution form are:

• -1 = normal error;

• 0 = log-normal error; and

• > 0 = Student’s t-distribution in natural log space with df equal to this value. For df > 30, results will be nearly identical to that for log-normal distribution. A df value of about 4 gives a fat-tail to the distribution. The se values entered in the data file must be the se in \(ln_{e}\) space.

Abundance indices typically assumed to have a log-normal error structure with units of se of \(ln_{e}\)(index). If the variance of the observations is available only as a cv (se of the observation divided by the mean value of the observation in natural space), then the value of standard error in natural log space can be calculated as \(\sqrt{(ln_e(1+(CV)^2))}\).

For the normal error structure, the entered values for se are interpreted directly as a se in arithmetic space and not as a cv. Thus switching from a log-normal to a normal error structure forces the user to provide different values for the se input in the data file.

If the data exist as a set of normalized Z-scores, you can assert a log-normal error structure after entering the data as \(exp(Z-score)\) because it will be logged by SS3. Preferably, the Z-scores would be entered directly and the normal error structure would be used.

7.10.0.3 Enable glssd Report

Indices with sd Report enabled will have the expected values for their historical values appear in the ss.std and ss.cor files. The default value is for this option is 0.

• 0 = sd Report not enabled for this index; and

• 1 = sd Report enabled for this index.

7.10.0.4 Data Format

• For fishing fleets, cpue is defined in terms of retained catch (biomass or numbers).

• For fishery independent surveys, retention/discard is not defined so cpue is implicitly in terms of total cpue.

• If a survey has its selectivity mirrored to that of a fishery, only the selectivity is mirrored so the expected cpue for this mirrored survey does not use the retention curve (if any) for the fishing fleet.

• If the fishery or survey has time-varying selectivity, then this changing selectivity will be taken into account when calculating expected values for the cpue or survey index.

• Year values that are before start year or after end year are excluded from model, so the easiest way to include provisional data in a data file is to put a negative sign on its year value.

• Duplicate survey observations for the same year are not allowed.

• Observations that are to be included in the model but not included in the negative log likelihood need to have a negative sign on their fleet ID. Previously the code for not using observations was to enter the observation itself as a negative value. However, that old approach prevented use of a Z-score environmental index as a “survey”. This approach is best for single or select years from an index rather than an approach to remove a whole index. Removing an index from the model should be done through the use of lambdas at the bottom of the control file which will eliminate the index from model fitting.

• Observations can be entered in any order, except if the super-year feature is used.

• Super-periods are turned on and then turned back off again by putting a negative sign on the season. Previously, super-periods were started and stopped by entering -9999 and the -9998 in the se field. See the Data Super-Period section of this manual for more information.

• If the statistical analysis used to create the cpue index of a fishery has been conducted in such a way that its inherent size/age selectivity differs from the size/age selectivity estimated from the fishery’s size and age composition, then you may want to enter the cpue as if it was a separate survey and with a selectivity that differs from the fishery’s estimated selectivity. The need for this split arises because the fishery size and age composition should be derived through a catch-weighted approach (to appropriately represent the removals by the fishery) and the cpue should be derived through an area-weighted approach to better serve as a survey of stock abundance.

7.11 Discard

If discard is not a feature of the model specification, then just a single input is needed:

If discard is being used, the input syntax is:

Note that although the user must specify a month for the observed discard data, the unit for discard data is in terms of a season rather than a specific month. So, if using a seasonal model, the input month values must correspond to some time during the correct season. The actual value will not matter because the discard amount is calculated for the entirety of the season. However, discard length or age observations will be treated by entered observation month.

7.11.0.1 Discard Units

The options are:

• 1 = values are amount of discard in either biomass or numbers according to the selection made for retained catch;

• 2 = values are fraction (in biomass or numbers) of total catch discarded, biomass/number selection matches that of retained catch; and

• 3 = values are in numbers (thousands) of fish discarded, even if retained catch has units of biomass.

7.11.0.2 Discard Error Distribution

The four options for discard error are:

• > 0 = df for Student’s t-distribution used to scale mean body weight deviations. Value of error in data file is interpreted as cv of the observation;

• 0 = normal distribution, value of error in data file is interpreted as cv of the observation;

• -1 = normal distribution, value of error in data file is interpreted as se of the observation;

• -2 = log-normal distribution, value of error in data file is interpreted as se of the observation in natural log space; and

• -3 = truncated normal distribution (new with v.3.30, needs further testing), value of error in data file is interpreted as se of the observation. This is a good option for low observed discard rates.

7.11.0.3 Discard Notes

• Year values that are before start year or after end year are excluded from model, so the easiest way to include provisional data in a data file is to put a negative sign on its year value.

• Negative value for fleet causes it to be included in the calculation of expected values, but excluded from the log likelihood.

• Zero (0.0) is a legitimate discard observation, unless log-normal error structure is used.

• Duplicate discard observations from a fleet for the same year are not allowed.

• Observations can be entered in any order, except if the super-period feature is used.

• Note that in the control file you will enter information for retention such that 1-retention is the amount discarded. All discard is assumed dead, unless you enter information for discard mortality. Retention and discard mortality can be either size-based or age-based (new with v.3.30).

7.11.0.4 Cautionary Note

The use of cv as the measure of variance can cause a small discard value to appear to be overly precise, even with the minimum se of the discard observation set to 0.001. In the control file, there is an option to add an extra amount of variance. This amount is added to the se, not to the cv, to help correct this problem of underestimated variance.

7.12 Mean Body Weight or Length

This is the overall mean body weight or length across all selected sizes and ages. This may be useful in situations where individual fish are not measured but mean weight is obtained by counting the number of fish in a specified sample (e.g., a 25 kg basket).

7.12.0.1 Partition

Mean weight data and composition data require specification of what group the sample originated from (e.g., discard, retained, discard + retained). Note: if retention is not defined in the selectivity section, observations with Partition = 2 will be changed to Partition = 0.

• 0 = combined catch in units of weight (whole, e.g., discard + retained);

• 1 = discarded catch in units of weight; and

• 2 = retained catch in units of weight.

7.12.0.2 Type

Specify the type of data:

• 1 = mean length; and

• 2 = mean body weight.

7.12.0.3 Observation - Units

Units must correspond to the units of body weight, normally in kg, (or mean length in cm). The expected value of mean body weight (or mean length) is calculated in a way that incorporates effect of selectivity and retention.

7.12.0.4 Error

Error is entered as the cv of the observed mean body weight (or mean length)

7.13 Population Length Bins

The first part of the length composition section sets up the bin structure for the population. These bins define the granularity of the alk and the coarseness of the length selectivity. Fine bins create smoother distributions, but a larger and slower running model. First read a single value to select one of three population length bin methods, then any conditional input for options 2 and 3:

7.13.0.1 Notes

There are some items for users to consider when setting up population length bins:

• For option 2, bin width should be a factor of min size and max size. For options 2 and 3, the data length bins must not be wider than the population length bins and the boundaries of the bins do not have to align. The transition matrix between population and data length bins is output to echoinput.sso.

• The mean size at settlement (virtual recruitment age) is set equal to the min size of the first population length bin.

• When using more, finer population length bins, the model will create smoother length selectivity curves and smoother length distributions in the alk, but run more slowly (more calculations to do).

• The mean weight-at-length, maturity-at-length and size-selectivity are based on the mid-length of the population bins. So these quantities will be rougher approximations if broad bins are defined.

• Provide a wide enough range of population size bins so that the mean body weight-at-age will be calculated correctly for the youngest and oldest fish. If the growth curve extends beyond the largest size bin, then these fish will be assigned a length equal to the mid-bin size for the purpose of calculating their body weight.

• While exploring the performance of models with finer bin structure, a potentially pathological situation has been identified. When the bin structure is coarse (note that some applications have used 10 cm bin widths for the largest fish), it is possible for a selectivity slope parameter or a retention parameter to become so steep that all the action occurs within the range of a single size bin. In this case, the model will see zero gradient of the log likelihood with respect to that parameter and convergence will be hampered.

• A value read near the end of the starter.ss file defines the degree of tail compression used for the alk, called alk tolerance. If this is set to 0.0, then no compression is used and all cells of the alk are processed, even though they may contain trivial (e.g., 1 e-13) fraction of the fish at a given age. With tail compression of, say 0.0001, the model, at the beginning of each phase, will calculate the min and max length bin to process for each age of each morph alk and compress accordingly. Depending on how many extra bins are outside this range, you may see speed increases near 10-20%. Large values of alk tolerance, say 0.1, will create a sharp end to each distribution and likely will impede convergence. It is recommended to start with a value of 0 and if model speed is an issue, explore values greater than 0 and evaluate the trade-off between model estimates and run time. The user is encouraged to explore this feature.

7.14 Length Composition Data Structure

If the value 0 is entered, then skip all length related inputs below and skip to the age data setup section. If value 1 is entered, all data weighting options for composition data apply equally to all partitions within a fleet. If the value 2 is entered, then the data weighting options are applied by the partition specified. Note that the partitions must be entered in numerical order within each fleet.

If the value for fleet is negative, then the vector of inputs is copied to all partitions (0 = combined, 1 = discard, and 2 = retained) for that fleet and all higher numbered fleets. This as a good practice so that the user controls the values used for all fleets.

7.14.0.1 Minimum Tail Compression

Compress tails of composition until observed proportion is greater than this value; negative value causes no compression; Advise using no compression if data are very sparse, and especially if the set-up is using age composition within length bins because of the sparseness of these data. A single fish being observed with tail compression on will cause the entire vector to be collapsed to that bin.

7.14.0.2 Added Constant to Proportions

Constant added to observed and expected proportions at length and age to make log likelihood calculations more robust. Tail compression occurs before adding this constant. Proportions are re-normalized to sum to 1.0 after constant is added.

The constant should be greater than 0. Commonly used values range from 0.00001 to 0.01. Larger values will cause differences among bins with smaller values to be less influential, leading to greater relative influence of the bins with the largest proportions of the compositions.

7.14.0.3 Combine Males & Females

Combine males into females at or below this bin number. This is useful if the sex determination of very small fish is doubtful so allows the small fish to be treated as combined sex. If Combine Males & Females > 0, then add males into females for bins 1 through this number, zero out the males, set male data to start at the first bin above this bin. Note that Combine Males & Females > 0 is entered as a bin index, not as the size associated with that bin. Comparable option is available for age composition data.

7.14.0.4 Compress Bins

This option allows for the compression of length or age bins beyond a specific length or age by each data source. As an example, a value of 5 in the compress bins column would condense the final five length bins for the specified data source.

7.14.0.5 Composition Error Distribution

The options are:

• 0 = Multinomial Error;

• 1 = Dirichlet-multinomial Error (linear); and

  • The Dirichlet-multinomial Error distribution requires the addition of a parameter lines for the natural log of the effective sample size multiplier (\(\theta\)) at the end of the selectivity parameter section in the control file. See the Dirichlet parameter in the control file for information regarding setup.

  • The Parameter Select option needs be used to specify which data sources should be weighted together or separate.

• 2 = Dirichlet-multinomial Error (saturation).

  • This parameterization of the Dirichlet-multinomial Error has not been tested, so this option should be used with caution. The Dirichlet-multinomial Error data weighting approach will calculate the effective sample size based on equation 12 from Thorson et al. (2017) where the estimated parameter will now be in terms of \(\beta\). The application of this method should follow the same steps detailed above for option 1.

7.14.0.6 Parameter Select

Value that indicates the groups of composition data for estimation of the Dirichlet parameter for weighting composition data.

• 0 = Default; and

• 1-N = Only used for the Dirichlet option. Set to a sequence of numbers from 1 to N where N is the total number of combinations of fleet and age/length. That is, if you have 3 fleets with length data, but only 2 also have age data, you would have values 1 to 3 in the length comp setup and 4 to 5 in the age comp setup. You can also have a data weight that is shared across fleets by repeating values in Parameter Select. Note that there can be no skipped numbers in the sequence from 1 to N, otherwise the model will exit on error when reading in the input files.

7.14.0.7 Minimum Sample Size

The minimum value (floor) for all sample sizes. This value must be at least 0.001. Conditional age-at-length data may have observations with sample sizes less than 1. Version 3.24 had an implicit minimum sample size value of 1.

7.14.0.8 Additional information on Dirichlet Parameter Number and Effective Sample Sizes

If the Dirichlet-multinomial error distribution is selected, indicate here which of a list of Dirichlet-multinomial parameters will be used for this fleet. So each fleet could use a unique Dirichlet-multinomial parameter, or all could share the same, or any combination of unique and shared. The requested number of Dirichlet-multinomial parameters are specified as parameter lines in the control file immediately after the selectivity parameter section. Please note that age-compositions Dirichlet-multinomial parameters are continued after length-compositions, so a model with one fleet and both data types would presumably require two new Dirichlet-multinomial parameters.

The Dirichlet estimates the effective sample size as \(N_{eff}=\frac{1}{1+\theta}+\frac{N\theta}{1+\theta}\) where \(\theta\) is the estimated parameter and \(N\) is the input sample size. Stock Synthesis estimates the natural log of the Dirichlet-multinomial parameter such that \(\hat{\theta}_{\text{fishery}} = e^{-0.6072} = 0.54\) where assuming \(N=100\) for the fishery would result in an effective sample size equal to 35.7.

This formula for effective sample size implies that, as the Stock Synthesis parameter \(ln(DM\text{\_theta})\) goes to large values (i.e., 20), then the adjusted sample size will converge to the input sample size. In this case, small changes in the value of the \(ln(DM\text{\_theta})\) parameter has no action, and the derivative of the negative log likelihood is zero with respect to the parameter, which means the Hessian will be singular and cannot be inverted. To avoid this non-invertible Hessian when the \(ln(DM\text{\_theta})\) parameter becomes large, turn it off while fixing it at the high value. This is equivalent to turning off down-weighting of fleets where evidence suggests that the input sample sizes are reasonable.

For additional information about the Dirichlet-multinomial please see Thorson et al. (2017) and the detailed Data Weighting section.

7.15 Length Composition Data

Composition data can be entered as proportions, numbers, or values of observations by length bin based on data expansions.

The data bins do not need to cover all observed lengths. The selection of data bin structure should be based on the observed distribution of lengths and the assumed growth curve. If growth asymptotes at larger lengths, having additional length bins across these sizes may not contribute information to the model and may slow model run time. Additionally, the lower length bin selection should be selected such that, depending on the size selection, to allow for information on smaller fish and possible patterns in recruitment. While set separately users should ensure that the length and age bins align. It is recommended to explore multiple configurations of length and age bins to determine the impact of this choice on model estimation.

Specify the length composition data as:

Note: the vector of length bins above will aggregate data from outside the range of values as follows:

Example of a single length composition observation:

7.15.0.1 Sex

If model has only one sex defined in the set-up, all observations must have sex set equal to 0 or 1 and the data vector by year will equal the number of the user defined data bins. This also applies to the age data.

In a 2 sex model, the data vector always has female data followed by male data, even if only one of the two sexes has data that will be used. The below description applies to a 2 sex model:

• Sex = 0 means combined male and female (must already be combined and information placed in the female portion of the data vector) (male entries must exist for correct data reading, then will be ignored).

• Sex = 1 means female only (male entries must exist for correct data reading, then will be ignored).

• Sex = 2 means male only (female entries must exist and will be ignored after being read).

• Sex = 3 means data from both sexes will be used and they are scaled so that they together sum to 1.0; i.e., sex ratio is preserved.

7.15.0.2 Partition

Partition indicates samples from either combined, discards, or retained catch. Note: if retention is not defined in the selectivity section, observations with Partition = 2 will be changed to Partition = 0.

• 0 = combined (whole, e.g., discard + retained);

• 1 = discard; and

• 2 = retained.

7.15.0.3 Excluding Data

• If the value of year is negative, then that observation is not transferred into the working array. This feature is the easiest way to include observations in a data file but not to use them in a particular model scenario.

• If the value of fleet in the length or age composition observed data line is negative, then the observation is processed and its expected value and log likelihood are calculated, but this log likelihood is not included in the total log likelihood. This feature allows the user to see the fit to a provisional observation without having that observation affect the model.

7.15.0.4 Note

When processing data to be input into SS3, all observed fish of sizes smaller than the first bin should be added to the first bin and all observed fish larger than the last bin should be condensed into the last bin.

The number of length composition data lines no longer needs to be specified in order to read the length (or age) composition data. Starting in v.3.30, the model will continue to read length composition data until a pre-specified exit line is read. The exit line is specified by entering -9999 at the end of the data matrix. The -9999 indicates to the model the end of length composition lines to be read.

Each observation can be stored as one row for ease of data management in a spreadsheet and for sorting of the observations. However, the 6 header values, the female vector and the male vector could each be on a separate line because admb reads values consecutively from the input file and will move to the next line as necessary to read additional values.

The composition observations can be in any order and replicate observations by a year for a fleet are allowed (unlike survey and discard data). However, if the super-period approach is used, then each super-periods’ observations must be contiguous in the data file.

7.16 Age Composition Option

The age composition section begins by reading the number of age bins. If the value 0 is entered for the number of age bins, then skips reading the bin structure and all reading of other age composition data inputs.

7.16.1 Age Composition Bins

If a positive number of age bins is read, then reads the bin definition next.

The bins are in terms of observed age (here age) and entered as the lower edge of each bin. Each ageing imprecision definition is used to create a matrix that translates true age structure into age structure. The first and last age’ bins work as accumulators. So in the example any age 0 fish that are caught would be assigned to the age = 1 bin.

7.16.2 Ageing Error

Here, the capability to create a distribution of age (e.g., age with possible bias and imprecision) from true age is created. One or many ageing error definitions can be created. For each, the model will expect an input vector of mean age and a vector of standard deviations associated with the mean age.

In principle, one could have year or laboratory specific matrices for ageing error. For each matrix, enter a vector with mean age for each true age; if there is no ageing bias, then set age equal to true age + 0.5. Alternatively, -1 value for mean age means to set it equal to true age plus 0.5. The addition of + 0.5 is needed so that fish will get assigned to the intended integer age. The length of the input vector is equal to the population maximum age plus one (0-max age), with the first entry being for age 0 fish and the last for fish of population maximum age even if the maximum age bin for the data is lower than the population maximum age. The following line is a vector with the standard deviation of age for each true age with a normal distribution assumption.

The model is able to create one ageing error matrix from parameters, rather than from an input vector. The range of conditions in which this new feature will perform well has not been evaluated, so it should be considered as a preliminary implementation and subject to modification. To invoke this option, for the selected ageing error vector, set the standard deviation of ageing error to a negative value for age 0. This will cause creation of an ageing error matrix from parameters and any age or size-at-age data that specify use of this age error pattern will use this matrix. Then in the control file, add a full parameter line below the cohort growth deviation parameter (or the movement parameter lines if used) in the mortality growth parameter section. These parameters are described in the control file section of this manual.

Code for ageing error calculation can be found in SS_miscfxn.tpl, search for function “get_age_age” or “SS_Label_Function 45”.

7.16.3 Age Composition Specification

If age data are included in the model, the following set-up is required, similar to the length data section. See Length Composition Data Structure for details on each of these inputs.

Syntax for Sex, Partition, and data vector are same as for length. The data vector has female values then male values, just as for the length composition data.

7.16.3.1 Age Error

Age error (Age Err) identifies which ageing error matrix to use to generate expected value for this observation.

7.16.3.2 Lbin Low and Lbin High

Lbin lo and Lbin hi are the range of length bins that this age composition observation refers to. Normally these are entered with a value of -1 and -1 to select the full size range. Whether these are entered as population bin number, length data bin number, or actual length is controlled by the value of the length bin range method above.

• Entering value of 0 or -1 for Lbin lo converts Lbin lo to 1;

• Entering value of 0 or -1 for Lbin hi converts Lbin hi to Maxbin;

• It is strongly advised to use the -1 codes to select the full size range. If you use explicit values, then the model could unintentionally exclude information from some size range if the population bin structure is changed.

• In reporting to the comp_report.sso, the reported Lbin_lo and Lbin_hi values are always converted to actual length.

7.16.3.3 Excluding Data

As with the length composition data, a negative year value causes the observation to not be read into the working matrix, a negative value for fleet causes the observation to be included in expected values calculation, but not in contribution to total log likelihood, a negative value for month causes start-stop of super-period.

7.17 Conditional Age-at-Length

Use of conditional age-at-length will greatly increase the total number of age composition observations and associated model run time, but there can be several advantages to inputting ages in this fashion. First, it avoids double use of fish for both age and size information because the age information is considered conditional on the length information. Second, it contains more detailed information about the relationship between size and age so provides stronger ability to estimate growth parameters, especially the variance of size-at-age. Lastly, where age data are collected in a length-stratified program, the conditional age-at-length approach can directly match the protocols of the sampling program.

However, simulation research has shown that the use of conditional age-at-length data can result in biased growth estimates in the presence of unaccounted for age-based movement when length-based selectivity is assumed (H. H. Lee et al. 2017), when other age-based processes (e.g., mortality) are not accounted for (H. Lee et al. 2019), or based on the age sampling protocol (Piner, Lee, and Maunder 2016). Understanding how data are collected (e.g., random, length-conditioned samples) and the biology of the stock is important when using conditional age-at-length data for a fleet.

In a two sex model, it is best to enter these conditional age-at-length data as single sex observations (sex = 1 for females and = 2 for males), rather than as joint sex observations (sex = 3). Inputting joint sex observations comes with a more rigid assumption about sex ratios within each length bin. Using separate vectors for each sex allows 100% of the expected composition to be fit to 100% observations within each sex, whereas with the sex = 3 option, you would have a bad fit if the sex ratio were out of balance with the model expectation, even if the observed proportion at age within each sex exactly matched the model expectation for that age. Additionally, inputting the conditional age-at-length data as single sex observations isolates the age composition data from any sex selectivity as well.

Conditional age-at-length data are entered within the age composition data section and can be mixed with marginal age observations for other fleets of other years within a fleet. To treat age data as conditional on length, Lbin_lo and Lbin_hi are used to select a subset of the total size range. This is different from setting Lbin_lo and Lbin_hi both to -1 to select the entire size range, which treats the data entered on this line within the age composition data section as marginal age composition data.

In this example observation, the age data is treated as on being conditional on the 2 cm length bins of 10–11.99, 12–13.99, 14–15.99, and 16–17.99 cm. If there are no observations of ages for a specific sex within a length bin for a specific year, that entry may be omitted.

7.18 Mean Length or Body Weight-at-Age

The model also accepts input of mean length-at-age or mean body weight-at-age. This is done in terms of observed age, not true age, to take into account the effects of ageing imprecision on expected mean size-at-age. If the value of the Age Error column is positive, then the observation is interpreted as mean length-at-age. If the value of the Age Error column is negative, then the observation is interpreted as mean body weight-at-age and the abs(Age Error) is used as Age Error.

7.18.0.1 Note

• Negatively valued mean size entries with be ignored in fitting. This feature allows the user to see the fit to a provisional observation without having that observation affect the model.

• A number of fish value of 0 will cause mean size value to be ignored in fitting. If the number of fish is zero, a non-zero mean size or body weight-at-age value, such as 0.01 or -999, still needs to be added. This feature allows the user to see the fit to a provisional observation without having that observation affect the model.

• Negative value for year causes observation to not be included in the working matrix. This feature is the easiest way to include observations in a data file but not to use them in a particular model scenario.

• Each sexes’ data vector and N fish vector has length equal to the number of age bins.

• The “Ignore” column is not used (set aside for future options) but still needs to have default values in that column (any value).

• Where age data are being entered as conditional age-at-length and growth parameters are being estimated, it may be useful to include a mean length-at-age vector with nil emphasis to provide another view on the model’s estimates.

• An experiment that may be of interest might be to take the body weight-at-age data and enter it to the model as empirical body weight-at-true age in the wtatage.ss file, and to contrast results to entering the same body weight-at-age data here and to attempt to estimate growth parameters, potentially time-varying, that match these body weight data.

• If using mean size-at-age data, please see the lambda usage notes regarding issues for model fitting depending upon other data within the model.

7.19 Environmental Data

The model accepts input of time series of environmental data. Parameters can be made to be time-varying by making them a function of one of these environmental time series. In v.3.30.16 the option to specify the centering of environmental data by either using the mean of the by mean and the z-score.

The final two lines in the example above indicate in that variable series 1 will be centered by subtracting the mean and dividing by the sd (indicated by the -1 value in the year column). The environmental variable series 2 will be centered by subtracting the mean of the time series (indicated by the -2 value in the year column). The input in the “value” column for both of the final two lines specifying the centering of the time series is ignored by the model. The control file also will need to be modified to in the long parameter line column “env-var” for the selected parameter. This feature was added in v.3.30.16.

7.19.0.1 Note

• Any years for which environmental data are not read are assigned a value of 0.0. None of the current link functions contain a link parameter that acts as an offset. Therefore, you should subtract the mean from your data. This lessens the problem with missing observations, but does not eliminate it. A better approach for dealing with missing observations is to use a different approach for the environmental effect on the parameter. Set up the parameter to have random deviations for all years, then enter the zero-centered environmental information as a special survey of type 35 and set up the catchability of that survey to be a link to the deviation vector. This is a more complex approach, but it is superior in treatment of missing values and superior in allowing for error in the environmental relationship.

• Users can assign environmental conditions for the initial equilibrium year by including environmental data for one year before the start year. However, this works only for recruitment parameters, not biology or selectivity parameters.

• Environmental data can be read for up to 100 years after the end year of the model. Then, if the recruitment-environment link has been activated, the future recruitments will be influenced by any future environmental data. This could be used to create a future “regime shift” by setting historical values of the relevant environmental variable equal to zero and future values equal to 1, in which case the magnitude of the regime shift would be dictated by the value of the environmental linkage parameter. Note that only future recruitment and growth can be modified by the environmental inputs; there are no options to allow environmentally-linked selectivity in the forecast years.

7.20 Generalized Size Composition Data

The generalized approach to size composition information was designed initially to provide a means to include weight frequency data. However, the uses are broader, such as allowing for size composition data with different data bins. The user can define as many generalized size composition methods as necessary.

• Each method has a specified number of bins.

• Each method has “units” so the frequencies can be in units of biomass or numbers.

• Each method has “scale” so the bins can be in terms of weight or length (including ability to convert bin definitions in pounds or inches to kg or cm).

• The composition data is input as females then males, just like all other composition data in SS3. In a two-sex model, the new composition data can be combined sex, single sex, or both sex.

• The generalized size composition data can be from the combined discard and retained (i.e., whole), discard only, or retained only.

• There are two options for treating fish that in population size bins are smaller than the smallest size frequency bin.

  • Option 1: By default, these fish are excluded (unlike length composition data where the small fish are automatically accumulated up into the first bin).

  • Option 2: If the first size bin is given as a negative value, then accumulation is turned on and the absolute value of the entered value is used as the lower edge of the first size bin.

Example input is shown below. Note that the format is identical to the length composition data, including sex and partition options, except for the addition of the first column, which indicates the size frequency method.

7.20.0.1 Note

• There is no tail compression for generalized size frequency data.

• Super-period capability is as for the length and age composition data.

• By choosing units = 2 and scale = 3 with identical bins and a negative first bin to turn accumulation of small fish on, the size composition method is identical to the length composition method.

• Bin boundaries do not need to align with the population length bin boundaries. The model interpolates as necessary.

• Size bins cannot be defined as narrower than the population bin width.

• The transition matrix can depend upon weight-at-length which differs between sexes and can vary seasonally. Thus, the transition matrix is calculated internally for each sex and each season.

7.21 Tag-Recapture Data

Each released tag group is characterized by an area, time, sex and age at release. Each recapture event is characterized by a time and fleet (since fleets operate in only one area, it is not necessary to specify the area of recapture). Fleets with tagging data must be fishing fleets (e.g., fleet type 1 or 2).

Inside the model, the tagged cohort is apportioned across all growth patterns in a given area at a given time (with options to apportion to only one sex or to both). The tag cohort by growth pattern then behaves according to the movement and mortality of the growth pattern. The number of tagged fish is modeled as a negligible fraction of the total population, so a tagging event does not move fish from an untagged group to a tagged group. Instead, tagged fish are seeded into the population with no impact at all on the total population abundance or mortality.

Predominant age at release for each tag group must be assigned; this requirement keeps SS3 efficient. By assigning a tag group to a single age rather than distributing it across all possible ages according to the size composition of the release group, the tag group can be tracked as a single cohort through the age by time matrix with minimal overhead to the rest of the model. Tags are released at the beginning of a season and recaptures follow the timing of the fleet that made the recapture.

7.21.0.1 Note

• The release data must be entered in tag group order.

• <tfill> values are placeholders and are replaced by program generated values for model time.

• Analysis of the tag-recapture data has one negative log likelihood component for the distribution of recaptures across areas and another negative log likelihood component for the decay of tag recaptures from a group over time. Note the decay of tag recaptures from a group over time suggests information about mortality is available in the tag-recapture data. More on this is in the control file documentation.

• Do tags option 2 adds an input compared to do tags option 1, minimum recaptures. Minimum recaptures option allows the user to exclude tag groups that have few recaptures after the mixing period from the likelihood. This may be useful when few tags from a group have been recaptured as an alternative to manually removing the groups with these low numbers of recaptured tags from the tagging data.

• Warning for earlier versions of SS3: A shortcoming in the recapture calculations when also using Pope’s \(F\) approach was identified and corrected in v.3.30.14.

7.22 Stock (Morph) Composition Data

It is sometimes possible to observe the fraction of a sample that is composed of fish from different stocks. These data could come from genetics, otolith microchemistry, tags, or other means. The growth pattern feature allows definition of cohorts of fish that have different biological characteristics and which are independently tracked as they move among areas. SS3 now incorporates the capability to calculate the expected proportion of a sample of fish that come from different growth patterns, “morphs”. In the inaugural application of this feature, there was a 3 area model with one stock spawning and recruiting in area 1, the other stock in area 3, then seasonally the stocks would move into area 2 where stock composition observations were collected, then they moved back to their natal area later in the year.

7.22.0.1 Note

• The number of stocks entered with these data must match the number of growth patterns (morphs) in the control file.

• Each data line for unique observations should enter data for morph 1 first followed sequentially for each morph included in the model.

• The expected value is combined across sexes. The entered data values will be normalized to sum to one within SS3.

• The “null” flag is included here in the data input section and is a reserved spot for future features.

• Note that there is a specific value of minimum compression to add to all values of observed and expected.

• Warning for earlier versions of SS3: A flaw was identified in the calculation of accumulation by morph. This has been corrected in v.3.30.14. Older versions were incorrectly calculating the catch by morph using the expectation around age-at-length which already was accounting for the accumulation by morph.

7.23 Selectivity Empirical Data (future feature)

It is sometimes possible to conduct field experiments or other studies to provide direct information about the selectivity of a particular length or age relative to the length or age that has peak selectivity, or to have a prior for selectivity that is more easily stated than a prior on a highly transformed selectivity parameter. This section provides a way to input data that would be compared to the specified derived value for selectivity. This is a placeholder at this time, required to include in the data file and will be fully implemented soon.

7.24 Excluding Data

Data that are before the model start year or greater than the retrospective year are not moved into the internal working arrays at all. So if you have any alternative observations that are used in some model runs and not in others, you can simply give them a negative year value rather than having to comment them out. The first output to data.ss_new has the unaltered and complete input data. Subsequent reports to data.ss_new produce expected values or bootstraps only for the data that are being used. Additional information on bootstrapping is available in Bootstrap Data Files Section.

Data that are to be included in the calculations of expected values, but excluded from the calculation of negative log likelihood, are flagged by use of a negative value for fleet number.

7.25 Data Super-Periods

The super-period capability allows the user to introduce data that represent a blend across a set of time steps and to cause the model to create an expected value for this observation that uses the same set of time steps. The option is available for all types of data and a similar syntax is used.

All super-period observations must be contiguous in the data file. All but one of the observations in the sequence will have a negative value for fleet ID so the data associated with these dummy observations will be ignored. The observed values must be combined outside the model and then inserted into the data file for the one observation with a positive fleet number.

Super-periods are started with a negative value for month, and then stopped with a negative value for month, observations within the super-period are designated with a negative fleet field. The standard error or input sample size field is now used for weighting of the expected values. An error message is generated if the super-period does not contain one observation with a positive fleet field.

An expected value for the observation will be computed for each selected time period within the super-period. The expected values are weighted according to the values entered in the se (or input sample size) field for all observations except the single observation holding the combined data. The expected value for that year gets a relative weight of 1.0. So in the example below, the relative weights are: 1982, 1.0 (fixed); 1983, 0.85; 1985, 0.4; 1986, 0.4. These weights are summed and rescaled to sum to 1.0, and are output in the echoinput.sso file.

Not all time steps within the extent of a super-period need be included. For example, in a three season model, a super-period could be set up to combine information from season 2 across 3 years, e.g., skip over the season 1 and season 3 for the purposes of calculating the expected value for the super-period. The key is to create a dummy observation (negative fleet value) for all time steps, except 1, that will be included in the super-period and to include one real observation (positive fleet value; which contains the real combined data from all the specified time steps).

A time step that is within the time extent of the super-period can still have its own separate observation. In the above example, the survey observation in 1984 could be entered as a separate observation, but it must not be entered inside the contiguous block of super-period observations. For composition data (which allow for replicate observations), a particular time steps’ observations could be entered as a member of a super-period and as a separate observation.

The super-period concept can also be used to combine seasons within a year with multiple seasons. This usage could be preferred if fish are growing rapidly within the year so their effective age selectivity is changing within year as they grow; fish are growing within the year so fishery data collected year round have a broader size-at-age modes than a mid-year model approximation can produce; and it could be useful in situations with very high fishing mortality.

8 Control File

8.1 Overview of Control File

These listed model features are denoted in the control file in the following order:

1. Number of growth patterns and platoons

2. Design matrix for assignment of recruitment to area/settlement event/growth pattern

3. Design matrix for movement between areas

4. Definition of time blocks that can be used for time-varying parameters

5. Controls far all time-varying parameters
   

6. Specification for growth and fecundity

7. Natural mortality growth parameters, weight-at-length, maturity, and fecundity, for each sex

8. Hermaphroditism parameter line (if used)

9. Recruitment distribution parameters for each area, settlement event, and growth pattern

10. Cohort growth deviation

11. Movement between areas (if used)

12. Age error parameter line (if used)

13. Catch multiplier (if used)

14. Fraction female

15. Setup for any mortality-growth parameters are time-varying

16. Seasonal effects on biology parameters
    

17. Spawner-recruitment parameters

18. Setup for any stock recruitment parameters are time-varying

19. Recruitment deviations
    

20. \(F\) ballpark value in specified year

21. Method for calculating fishing mortality (\(F\))

22. Initial equilibrium \(F\) for each fleet
    

23. Catchability (Q) setup for each fleet and survey

24. Catchability parameters

25. Setup for any catchability parameters are time-varying
    

26. Length selectivity, retention, discard mortality setup for each fleet and survey

27. Age selectivity setup for each fleet and survey

28. Parameters for length selectivity, retention, discard mortality for each fleet and survey

29. Parameters for age selectivity, retention, discard mortality for each fleet and survey

30. Setup for any selectivity parameters that are time-varying
    

31. Tag-recapture parameters
    

32. Variance adjustments

33. Lambdas for likelihood components

The order in which they appear in the control file has grown over time rather opportunistically, so it may not appear particularly logical at this time, especially various aspects of recruitment distribution and growth.

8.2 Parameter Line Elements

The primary role of the control file is to define the parameters to be used by the model. The general syntax of the 14 elements of a long parameter line is described here. If used, time-varying parameter lines use only the first seven elements of a parameter line and will be referred to as a short parameter line. Three types of time-varying properties can be applied to a base parameter: blocks or trend, environmental linkage, and random deviation. Each parameter line contains:

Note that relative to Stock Synthesis v.3.24, the order of PRIOR sd and PRIOR TYPE have been switched and the PRIOR TYPE options have been renumbered.

The full parameter line (14 in length) syntax for the mortality-growth, spawn-recruitment, catchability, and selectivity sections provides additional controls to give the parameter time-varying properties. If a parameter (a full parameter line of length 14) is set up to be time-varying (i.e., parameter time blocks, annual deviations), short parameter lines, the first 7 elements, are required to be specified immediately after the main parameter block (i.e., mortality-growth parameter section). Additional information regard time-varying parameters and how to implement them is in the Using Time-Varying Parameters section.

8.3 Terminology

The term COND appears in the “Typical Value” column of this documentation (it does not actually appear in the model files), it indicates that the following section is omitted except under certain conditions, or that the factors included in the following section depend upon certain conditions. In most cases, the description in the definition column is the same as the label output to the ss_new files.

8.4 Beginning of Control File Inputs

8.4.1 Weight-at-Age

The capability to read empirical body weight-at-age for the population and each fleet was added starting in v.3.04, in lieu of generating these weights internally from the growth parameters, weight-at-length, and size-selectivity. The values are read from a separate file named, wtatage.ss. This file is only required to exist if this option is selected. See the Empirical Weight-at-Age section for additional information on file formatting for empirical weight-at-age.

8.4.2 Settlement Timing for Recruits and Distribution

In older versions of SS3 one value of spawning biomass was calculated annually at the beginning of one specified spawning season and this spawning biomass produced one annual total recruitment value. The annual recruitment value was then distributed among seasons, areas, and growth types according to other model parameters.

Additional control of the seasonal timing was added in v.3.30 and now there is an explicit elapsed time between spawning and recruitment. Spawning still occurs, just once per year, which defines a single spawning biomass for the stock-recruitment curve, but its timing can be at any specified time, not just the beginning of a season. Recruitment of the progeny from an annual spawning can now enter the population in one or more settlement events, at some point after spawning as defined by the user.

The above example specifies settlement to mid-May (month 5.5). Note that normally the calendar age at settlement is 0 if settlement happens between the time of spawning and the end of that year, and at age 1 if settlement is in the year after spawning.

Below is an example setup where there are multiple settlement events, with one occurring the following year after spawning:

Details regarding settlement of recruits and timing:

• Recruitment happens in specified settlement events (growth pattern, month, area).

• Number of unique settlement timings is calculated at runtime.

• Now there is explicit elapsed time between spawning and recruitment.

• Growth and natural mortality of the platoon begins at time of settlement, which is its real age 0.0 for growth; but pre-settlement fish exist from the beginning of the season of settlement, so can be caught if selected.

• Age at recruitment now user-controlled (should be 0 if in year of spawning).

• All fish become integer age 1 (for age determination) on their first January 1st.

• Recruitment can occur > 12 months after spawning which is achieved by setting the settlement age to a value greater than 1.0.

The distribution of recruitment among these settlement events is controlled by recruitment apportionment parameters. There must be a parameter line for each growth pattern, then for each area, then for each settlement. All of these are required, but only those growth pattern \(\times\) area \(\times\) settlements designated to receive recruits in the recruitment design matrix will have the parameter used in the recruitment distribution calculation. For the recruitment apportionment, the parameter values are the natural log of apportionment weight. The sum of all apportionment weights is calculated for each growth pattern \(\times\) area \(\times\) settlements that have been designated to receive recruits in the recruitment design matrix. Then the apportionment weights are scaled to sum to 1.0 so that the total recruitment from the spawning event is distributed among the cells designated to receive recruitment. Additionally, these distribution parameters can be time-varying, so the fraction of the recruits that occur in a particular growth pattern, area, or settlement can change from year to year. To specify annual variation in the distribution or recruits by area add a start and end year in the deviation min year and max year columns. Similar to the apportionment of recruits by area, one should be fixed while the other area(s) can deviate relative to the one area. If annual deviations are specified then two additional short parameter lines will be required to specify the standard error and the autocorrelation for each area with deviations.

8.4.2.1 Recruitment Distribution and Parameters

Recruits are apportioned according to:

\[\text{apportionment}_i = \frac{e^{p_i}}{\sum_{i=1}^{N}e^{p_i}}\]

where \(p_i\) is the proportion of recruits to area \(i\) and \(N\) is the number of settlement events. These parameters are defined in the mortality-growth parameter section.

Tips for fixing or estimating the recruitment apportionment:

• Set the value for one of these parameters, \(p_i\), to 0.0 and not estimate it so that other parameters will be estimated (if not fixed) relative to its fixed value.

• Give the estimated parameters a min-max so they have a good range relative to the base parameter (i.e., of min = -5 and max = 5).

• In order to get a different distribution of recruitments in different years, you will need to make at least one of the recruitment distribution parameters time-varying.

In a seasonal model, all cohorts graduate to the age of 1 when they first reach January 1, even if the seasonal structure of the model has them being spawned in the late fall. In general, this means that the model operates under the assumption that all age data have been adjusted so that fish are age 0 at the time of spawning and all fish graduate to the next age on January 1. This can be problematic if the ageing structures deposit a ring at another time of year. Consequently, you may need to add or subtract a year to some of your age data to make it conform to the model expected data structure, or more ideally you may need to define the calendar year within the model to start at the beginning of the season at which ring deposition occurs. Talk with your ageing lab about their criteria for seasonal ring deposition.

Seasonal recruitment is coded to work smoothly with growth. If the recruitment occurring in each season is assigned the same growth pattern, then each seasonal cohort’s growth trajectory is simply shifted along the age/time axis. At the end of the year, the early born cohorts will be larger, but all are growing with the same growth parameters, so all will converge in size as they approach their common maximum length (e.g., no seasonal effects on growth).

At the time of settlement, fish are assigned a size equal to the lower edge of the first population size bin, and they grow linearly until they reach the age A1. A warning is generated if the first population length bin is greater than 10 cm as this seems an unreasonably large value for a larval fish. A1 is in terms of real age elapsed since birth. All fish advance to the next integer age on January 1, regardless of birth season. For example, consider a 2 season model with some recruitment in each season and with each season’s recruits coming from the same GP. At the end of the first year, the early born fish will be larger but both of the seasonal cohorts will advance to an integer age of 1 on Jan 1 of the next year. The full growth curve is still calculated below A1, but the size-at-age used is the linear replacement. Because the linear growth trajectory can never go negative, there is no need for the additive constant to the standard deviation (necessary for the growth model used in SS2 V1.x), but the option to add a constant has been retained in the model.

8.4.3 Movement

Here the movement of fish between areas are defined. This is a box transfer with no explicit adjacency of areas, so fish can move from any area to any other area in each time step. While not incorporated yet, there is a desire for future versions of SS3 to have the capability to allow sex-specific movement, and also to allow some sort of mirroring so that sexes and growth patterns can share the same movement parameters if desired.

Two parameters will be entered later for each growth pattern, area pair, and season.

• Movement is constant at the first parameter (P1) below the specified minimum age for movement change, constant at the second parameter (P2) above maximum age for movement change, and linearly interpolated for intermediate ages.

• A movement rate parameter can be set to use the same value as the corresponding parameter for the first defined movement pattern by entering a parameter value of -9999 and a negative phase value.

• For each source area, the implicit movement parameter value is 0.0 (movement within a single area). However, this default value can be replaced if the stay movement is selected to have an explicit pair of parameters (e.g., specify movement rate for area 1 to area 1) and will require additional parameter lines.

• A constant movement rate across all ages can be accomplished by either:

  • Setting both movement ages to 0, not estimating the first movement parameter, and using a second movement parameter to cover all ages from 0 to the maximum number of ages.

  • Setting movement ages to any value, estimating the first movement parameter, and setting the second movement parameter to have a value of -9998 with a negative phase.

• The parameter is exponentiated so that a movement parameter value of 0 becomes 1.0.

• For each source area, all movement rates are then summed and divided by this sum so that 100% of the fish are accounted for in the movement calculations. \[\text{rate}_i = \frac{e^{p_i}}{\sum_{j=1}^{N}e^{p_i}}\]

• At least one movement parameter must be fixed so that all other movement parameters are estimated relative to it. This is achieved naturally by not specifying the stay rate parameter, so it has a fixed value of 0.0.

• The resultant movement rates are multiplied by season duration in a seasonal model.

8.4.4 Time Blocks

Blocks and other time-vary parameter controls are operative during forecast years, so care should be taken when setting the end year of the last block in a pattern. If that end year is set to the last year in the time series, then the parameter will revert to the base value for the forecast. If the user wants to continue the last block through the forecast, it is advisable to set the last block’s end year value to -2 to cause SS3 to reset it to the last year of the forecast. Using the value -1 will set the block’s end year to the last year of the time series and leave the forecast at the base parameter value. Note that additional controls on time-varying parameters in forecast years are in the forecast section.

8.4.5 Auto-generation

Auto-generation is a useful way to automatically create the required short time-varying parameter lines which will be written in the control.ss_new file. These parameter lines can then be copied into the control file and modified as needed. As example, if you want to add a block to natural mortality, modify the block and block function entry of the mortality parameter line, ensure that auto-generation is set to 0 (for the biology section at least) and run the model without estimation. The control.ss_new file will now show the required block parameter line specification for natural mortality and this line can be copied into the main control file. Note, that if auto-generation is on (set to 0), the model will not expect to read the time-varying parameters in that section of the control file and will error out if they are present

8.5 Biology

8.5.1 Natural Mortality

Natural mortality (\(M\)) options include some options that are referenced to integer age and other options to real age since settlement. If using an option that references \(M\) to real age since settlement, \(M\) varies by age and will change by season (e.g., cohorts born early in the year will have different \(M\) than cohorts born later in the year).

8.5.1.1 Lorenzen Natural Mortality

Lorenzen natural mortality is based on the concept that natural mortality is driven by physiological and ecological processes and varies over the life cycle of a fish. So, natural mortality is scaled by the length of the fish. In this implementation, a reference age and \(M\) value are read in, and other ages will have an \(M\) scaled to its body size-at-age. However, if platoons are used, all will have the same \(M\) as their growth pattern. Lorenzen \(M\) calculation will be updated if the starting year growth parameters are active, but if growth parameters vary during the time series, the \(M\) is not further updated. Additionally, the \(M\) is linked to the length-at-age from the growth parameters and can’t be used for an empirical weight-at-age model. Be careful in using Lorenzen when there is time-varying growth.

8.5.1.2 Age-specific \(M\) Linked to Age-Specific Length and Maturity

This is an experimental option available as of v.3.30.17.

A general model for age- and sex-specific natural mortality expands a model developed by Mark N. Maunder et al. (2010) and Mark N. Maunder (2011) and is based on the following some assumptions:

1. \(M\) for younger fish is due mainly to processes that are functions of the size of the individuals (e.g., predation);

2. \(M\) increases after individuals become reproductively mature;

3. Maturity follows a logistic curve; and

4. \(M\) caused by senescence is either small or occurs at an age for which there are few fish alive, so it is not influential.

The model is based on combining the observation that \(M\) is inversely proportional to length for young fish (Kai Lorenzen 2000) and the logistic model from Lehodey, Senina, and Murtugudde (2008) for older fish. Natural mortality for a given sex and age is: \[M_{a,s} = M_{juv,s}\frac{L_{a,s}}{L_{mat*,s}}^{\lambda} + \frac{M_{mat,s}-M_{juv,s}\frac{L_{a,s}}{L_{mat*,s}}^{\lambda}}{1+e^{\beta_s(L_{a,s}- L_{50,s})}},\]

where \(M_{juv,s}\) (juvenile natural mortality), \(\lambda\) (power), \(L_{mat*,s}\) (first mature length of fish), and \(M_{mat,s}\) (the mature instantaneous natural mortality rate by sex, are user inputs in long parameter lines. For sub-option 1, \(L_{50}\) and \(\beta\) (slope) parameters taken from the maturity relationship within the model, which must use maturity-fecundity option 1. For sub-option 3, the \(L_{50,s}\) (the length at which 50% of fish are mature) and \(\beta\) (slope) parameters are specified in long parameter lines by the user.

Note that juvenile natural mortality, \(M_{juv,s}\), and first mature length of fish, \(L_{mat*,s}\), inputs are by sex (and growth pattern), but it is recommended to share them across sex by using the offset option. Using offset option 2 (males offset from females) causes male parameters to be an offset to the female parameters, so a parameter value of 0.0 for a male parameter will fix the parameter as same as the female parameter. Alternatively, using offset option 1 and setting males to 0.0 and not estimating the parameter fixes the parameter at the value of the female parameter (the section on fixing male parameters the same as female parameters has more details). This fulfills an additional assumption: \(M\) caused by reproduction may differ by sex, but juvenile \(M\) is independent of sex.

The length for a given age and sex, \(L_{a,s}\) is calculated within the model.

Some suggested defaults for user-provided parameter inputs are:

• \(\lambda = -1.5\) from Gulland (1987)

• \(M_{mat,s}=\frac{5.4}{t_{max,s}}\) from Hamel (submitted) if \(t_{max}\) is available, otherwise \(M_{mat,s} = 4.118K_{s}^{0.73}L_{inf,s}^{-0.33}\) as in Then et al. (2015)

• \(M_{juv,s} = 3W_{mat}^{-0.288}\) from K. Lorenzen (1996)

8.5.1.3 Age-range Lorenzen

The original implementation of Lorenzen natural mortality in Stock Synthesis uses a reference age and its associated natural mortality as inputs to determine the Lorenzen curve. However, sometimes this information is not known. The age-range Lorenzen instead uses a range of ages and the average natural mortality over them to calculate a Lorenzen natural mortality curve.

Like the original Lorenzen options, ages will have an \(M\) scaled to its body size-at-age and care should be taken when there are multiple growth patterns or time-varying growth.

8.5.1.4 Natural Mortality Options

8.5.2 Growth

8.5.2.1 Timing

When fish recruit at the real age of 0.0 at settlement, they have body size equal to the lower edge of the first population size bin. The fish then grow linearly until they reach a real age equal to the input value “growth-at-age for L1” and have a size equal to the parameter value for L1 (the minimum length parameter). As they age further, they grow according the selected growth equation. The growth curve is calibrated to go through the size L2 parameter when they reach the age of maximum length.

8.5.2.2 Maximum Length (Linf)

If “Growth at age for L2” is set equal to 999, then the size at the L2 parameter is used as Linf.

8.5.2.3 von Bertalanffy growth function

The von Bertalanffy growth curve is parameterized as:

\[L_t = L_\infty + (L_{1}-L_\infty)e^{-k(a-A_{1})}\]

with parameters \(L_{1}\), \(L_\infty\), and \(k\). The \(L_\infty\) is calculated as:

\[L_\infty = L_{1} + \frac{(L_2 - L_1)}{1-e^{-k(A2-A1)}}\]

based on the input values of fixed age for first size-at-age (\(A_1\)) and fixed age for second size-at-age (\(A_2\)).

8.5.2.4 Richards growth function

The Richards (1959) growth model as parameterized by Schnute (1981) provides a flexible growth parameterization that allows for a variety of growth curve shapes. The Richards growth is invoked by entering option 2 in the growth type field. The Richards growth function uses the standard growth parameters (\(L_1\), \(L_2\), \(k\)) and a fourth shape parameter \(b\) that is specified after the growth coefficient \(k\).

The Richards growth model is parameterized as:

\[L_t = \left[L_1^b + (L_2^b-L_1^b)\frac{1-e^{-k(t-A_{1})}}{1-e^{-k(A_2-A_1)}}\right]^{1/b}\]

with parameters \(L_1\), \(L_2\), \(k\), and \(b\).

The \(b\) shape parameter can be positive or negative but not precisely 0. When estimating \(b\) as a floating-point number, there is effectively no risk of the parameter becoming precisely zero during estimation, as long as the initial value is non-zero.

As special cases of the Richards growth model, \(b\!=\!1\) is von Bertalanffy growth and \(b\) near 0 is Gompertz growth. To use a Gompertz growth curve, the \(b\) parameter can be fixed at a small value such as 0.0001.

When \(A_1\) is greater than the youngest age in the model, some combinations of Richards growth parameters can lead to undefined (NaN) predicted length for the younger ages. The choice of \(A_1\) and \(A_2\) will affect the possible growth curve shapes.

The SS3 website includes a vignette providing further technical insights for using the Richards growth model in Stock Synthesis.

8.5.2.5 Mean size-at-maximum age

The mean size of fish in the max age bin depends upon how close the growth curve is to Linf by the time it reaches max age and the mortality rate of fish after they reach max age. Users specify the mortality rate to use in this calculation during the initial equilibrium year. This must be specified by the user and should be reasonably close to \(M\) plus initial \(F\). In v.3.30, this uses the von Bertalanffy growth out to 3 times the maximum population age and decays the numbers at age by exp(-value set here). For subsequent years of the time series, the model should update the size-at-maximum age according to the weighted average mean size of fish already at maximum age and the size of fish just graduating into maximum age. Unfortunately, this updating is only happening in years with time-varying growth. This will hopefully be fixed in the future version.

8.5.2.6 Age-specific K

This option creates age-specific K multipliers for each age of a user-specified age range, with independent multiplicative factors for each age in the range and for each growth pattern/sex. The null value is 1.0 and each age’s K is set to the next earlier age’s K times the value of the current age’s multiplier. Each of these multipliers is entered as a full parameter line, so inherits all time-varying capabilities of full parameters. The lower end of this age range cannot extend younger than the specified age for which the first growth parameter applies. This is a beta model feature, so examine output closely to assure you are getting the size-at-age pattern you expect. Beware of using this option in a model with seasons within year because the K deviations are indexed solely by integer age according to birth year. There is no offset for birth season timing effects, nor is there any seasonal interpolation of the age-varying K.

8.5.2.7 Growth cessation

A growth cessation model was developed for the application to tropical tuna species (Mark N. Maunder et al. 2018). Growth cessation allows for a linear relationship between length and age, followed by a marked reduction of growth after the onset of sexual maturity by assuming linear growth for the youngest individuals and then a logistic function to model the decreasing growth rate at older ages.

8.5.3 Maturity-Fecundity

8.5.4 Hermaphroditism

Sequential hermaphroditism can be modeled in Stock Synthesis by having recruits be all female or all male and then defining a 3-parameter function which represents the annual age-specific probability of transition to the other sex. In a seasonal model, the transition can occur after each season or just once per year. There are also settings to control the first age which transitions and the contribution of males to the spawning biomass calculations.

The fraction female parameter should be configured to reflect the fraction of age 0 fish of each sex. Values of 0 or 1 may lead to NaN likelihoods, so we recommend changing values that would be fixed at 0 or 1 to 0.000001 or 0.999999.

The hermaphroditism option requires three full parameter lines in the mortality and growth section. These parameters control a cumulative normal distribution function as follows:

1. The inflection age where the transition rate is halfway from 0 to its asymptote.

2. The standard deviation of the cumulative normal (such that about 95% of the increase in transition rate occurs within 2 standard deviations of the inflection point).

3. The asymptotic transition rate (the highest proportion that will transition to the other sex in a time step).

These parameter lines are entered directly after the weight-at-length parameters for males.

8.5.5 Natural Mortality and Growth Parameter Offset Method

The most common setup for natural mortality and growth parameters for two-sex models is direct assignment (option 1) which allows for independent estimation (or fixing) of natural mortality and growth parameters by sex. Within the direct assignment option there is functionality to set male parameters equal to the corresponding female parameter if the male INIT value is set to 0 and the phase is negative.

Alternatively, there may be situations where a user wants to create direct linkages between natural mortality and growth parameters between sexes (options 2 or 3). If the parameter offset option 2 is selected, the control file still requires that all male natural mortality and growth parameters lines to be included. The natural mortality and growth parameters (e.g., k, Lmin, Lmax, CV1, CV2) for sex > 1 (typically male fish) have a value that is an exponential offset to the female natural mortality and growth parameters, e.g., \(M_{\text{male}} = M_{\text{female}}*exp(M_{\text{male offset}})\). An offset parameter can be fixed at 0.0, at a non-zero value, or estimated.

Parameter offset option 3 has an offset feature for the growth cv and for natural mortality. For the growth cv, the parameter for cv at old age is an exponential offset from the parameter for cv at young age, e.g., \(CV_{\text{old}} = CV_{\text{young}}*exp(CV_{\text{offset}})\). This allows for cv old to track an estimated cv young parameter. For natural mortality, if there is more than 1 natural mortality parameter, then parameters 2 and higher for the same sex and growth pattern are exponential offsets from the first natural mortality parameter. Note that it is an old feature designed to work with natural mortality option 1 (breakpoints). It may work with natural mortality options 3 and 4, but this has not been tested.

8.5.6 Catch Multiplier

These parameter lines are only included in the control file if the catch multiplier field in the data file is set to 1 for a fleet. The model expected catch \(C_{exp}\) by fleet is estimated by:

\[C_{exp} = \frac{C_{obs}}{c_{mult}}\]

where \(C_{obs}\) is the input catch by fleet (observed catch) within the data file and \(c_{mult}\) is the estimated (or fixed) catch multiplier. It has year-specific, not season-specific, time-varying capabilities. In the catch likelihood calculation, expected catch is multiplied by the catch multiplier by year and fishery to get \(C_{obs}\) before being compared to the observed retained catch as modified by the \(c_{mult}\).

8.5.7 Ageing Error Parameters

These parameters are only included in the control file if one of the ageing error definitions in the data file has requested this feature (by putting a negative value for the ageing error of the age zero fish of one ageing error definition). As of v.3.30.12, these parameters now have time-varying capability. Seven additional full parameter lines are required. The parameter lines specify:

1. Age at which the estimated pattern begins (just linear below this age), this is the start age.

2. Bias at start age (as additive offset from unbiased age).

3. Bias at maximum (as additive offset from unbiased age).

4. Power function coefficient for interpolating between those 2 values (value of 0.0 produces linear interpolation in the bias).

5. Standard deviation at start age.

6. Standard deviation at max age.

7. Power function coefficient for interpolating between those 2 values.

Code for implementing vectors of mean age and standard deviation of age can be located online within the SS_miscfxn.tpl file, search for function “get_age_age” or “SS_Label_Function 45”.

8.5.8 Sex Ratio

The last line in the mortality-growth parameter section allows the user to fix or estimate the sex ratio between female and male fish at recruitment. The parameter is specified in the fraction of female fish and is applied at settlement. The default option is a sex ratio of 0.50 with this parameter not being estimated. Any composition data input as type = 3, both sexes, will be informative to the sex ratio because it scales females and males together, not separately, for this data type input. Estimation of the sex ratio is a new feature and should be done with care with the user checking that the answer is reflective of the data.

As of v.3.30.12, this parameter now has time-varying capability similar to other parameters in the mortality-growth section.

8.5.9 Predator Fleet Mortality

The ability to define a predator fleet was first implemented in v.3.30.18. A parameter line for predator mortality is only required if a predator fleet has been defined in the data file. For each fleet that is designated as a predator, a new parameter line is created in the mg section in the control file. This parameter will have the label \(M2\text{\_pred1}\), where the “1” is the index for the predator (not the index of the fleet being used as a predator). More than one predator can be included. If the model has > 1 season, it is normal to expect \(M2\) to vary seasonally. Therefore, only if the number of seasons is greater than 1, follow each \(M2\) parameter with number of season parameters to provide the seasonal multipliers. These are simple multipliers times \(M2\), so at least one of these needs to have a non-estimated value. The set of multipliers can be used to set \(M2\) to only operate in one season if desired. If there is more than one predator fleet, each will have its own seasonal multipliers. If there is only 1 season in the model, then no multiplier lines are included.

\(M2\) is age-specific, but not sex or morph specific. The value of the \(M2\) parameter will be distributed across ages according to the selectivity for this fleet. In this example note that “pred1” refers to the first predator in the model, note the fleet number in which that predator has been 3 configured. The resultant age-specific \(M2\) is added to the base \(M\) to create a total age-specific \(M\) that operates in the model exactly as \(M\) has always operated.

Because \(M2\) is a mg, it can be time-varying like any other mg. This is important because \(M2\), as a component added to base \(M\), will probably always need to be time-varying by blocks, random walk or linkage to external driver. A time series of \(M2\) from an external source could be input by setting the \(M2\) parameter to have a base value of 0.0 and linking to the time series in the environmental data section of the data file using an additive link. In addition, the relationship should have a fixed slope of 1.0 such that \(M2(y) = 0.0 + 1.0 * M2\text{\_env\_input}(y)\).

Note that all existing reports of natural mortality are the total (base \(M\) + \(M2\)) natural mortality. The \(M2\) parameter is active in the virgin year and initial equilibrium year, where the value of \(M2\) in the start year is used. In the future, separate control of \(M2\) for the initial equilibrium will be provided. \(M2\) is part of the total \(M\) used in the spr and msy benchmark calculations. \(M2\) is active in the forecast era, so be attentive to its configuration if it is time-varying. Testing to date shows that this \(M2\) feature can replicate previous results using bycatch fleets.

Note that predator calculations probably will fail if tried with \(F\) Method = 1 (Pope’s), even though the Pope calculation is built into the hybrid \(F\) approach.

8.5.10 Read Biology Parameters

Next, the model reads the mg parameters in generally the following order (may vary based on selected options):

Example format for mortality-growth parameter section with 2 sexes, 2 areas. Parameters marked with COND are conditional on selecting that feature:

8.5.10.1 Setting Male Parameters Equal to Females

The model allows a short-cut for males to use the same parameter values as female fish for natural mortality, length minimum (Length at Amin), maximum length (Length at Amax), coefficient at younger ages (CV1), coefficient at older ages (CV2), and the growth coefficient (K) when using offset option = 1 (the offset section has information on the options available). If the INIT parameter value for males is set equal to 0.0 and the phase set to negative, not estimated, each of these male parameters will use the corresponding female parameter value for the males.

8.5.11 Time-varying Parameters

Please see the Time-Varying Parameter Specification and Setup section for details on how to set up time varying parameters. In short, additional short parameter lines will be needed after the long parameter lines. There are some additional considerations for time-varying growth.

8.5.12 Seasonal Biology Parameters

Seasonal effects are available for weight-length parameters, maturity, fecundity, and for the growth parameter K. The seasonal parameter values adjust the base parameter value for that season. \[P'=P*exp(\text{seas\_value})\]

8.6 Spawner-Recruitment

The spawner-recruitment section starts by specification of the functional relationship that will be used.

8.6.0.1 Equilibrium Recruitment

In principle, steepness should always be used when calculating equilibrium recruitment. This was not the default in early versions of Stock Synthesis, so has not come into common practice. The original logic, from early 1990s version of Stock Synthesis for long-lived U.S. west coast groundfish, was that fishing had not yet gone on long enough to have reduced spawning biomass enough to reduce expected recruitment noticeably for the chosen initial equilibrium year.

Steepness should be used in the equilibrium calculation whenever you believe that the initial equilibrium catch was large enough to have reduced expected recruitment below \(R_{0}\). Note that when using this option, the initial equilibrium catch cannot be greater than msy. SS3 uses the identical code for initial equilibrium catch and for msy calculation, so it is axiomatic that msy will be \(<=\) initial equilibrium catch. So, SS3 will estimate a \(R_{0}\) large enough to make msy < initial equilibrium catch. Alternatively, you can elect to not use this option and instead add many years to the beginning of the time series with that same level of initial catch. In this case, it is not equilibrium, so the catch can reduce the population below \(B_{MSY}\).

8.6.1 Spawner-Recruitment Functions

The number of age-0 fish is related to spawning biomass according to a stock-recruitment relationship. There are a number of options for the shape of the spawner-recruitment relationship: Beverton-Holt, Ricker, Hockey-Stick, and a survival-based stock recruitment relationship.

8.6.1.1 Beverton-Holt

The Beverton-Holt Stock Recruitment curve is calculated as: \[{R_y = \frac{4hR_0SB_y}{SB_0(1-h)+SB_y(5h-1)}e^{-0.5b_y\sigma^2_R+\tilde{R}_y}\qquad \tilde{R}_y\sim N(0;\sigma^2_R)}\]

where \(R_0\) is the unfished equilibrium recruitment, \(SB_0\) is the unfished equilibrium spawning biomass (corresponding to \(R_0\)), \(SB_y\) is the spawning biomass at the start of the spawning season during year \(y\), \(h\) is the steepness parameter, \(b_y\) is the bias adjustment fraction applied during year \(y\), is the standard deviation among recruitment deviations in natural log space, and is the log-normal recruitment deviation for year \(y\). The bias-adjustment factor (Methot and Taylor 2011) ensures unbiased estimation of mean recruitment even during data-poor eras in which the maximum likelihood estimate of the recruitment deviation is near 0.0.

8.6.1.2 Ricker

The Ricker Stock Recruitment curve is calculated as: \[{R_y = \frac{R_0SB_y}{SB_0}e^{h(1-SB_y/SB_0)}e^{-0.5b_y\sigma^2_R+\tilde{R}_y}\qquad \tilde{R}_y\sim N(0;\sigma^2_R)}\]

where the stock recruitment parameters have the same meaning as described above for the Beverton-Holt.

8.6.1.3 Hockey-Stick

The hockey-stick recruitment curve is calculated as: \[{R_y = join(R_{\text{min}}R_0+R_0\frac{SB_y}{hSB_0}(1-R_{\text{min}}))+R_0(1-join)e^{-0.5b_y\sigma^2_R+\tilde{R}_y}\qquad \tilde{R}_y\sim N(0;\sigma^2_R)}\] where \(R_{\text{min}}\) is the minimum recruitment level predicted at a spawning size of zero and is set by the user in the control file, \(h\) is defined as the fraction of \(SB_0\) below which recruitment declines linearly, and \(join\) is defined as: \[{ join = \bigg[1+e^{1000*\frac{(SB_0-hSB_0)}{SB_0}}\bigg]^{-1} }\]

8.6.1.4 Survivorship

The survivorship stock recruitment relationship based on Taylor et al. (2013) is a stock-recruitment model that enables explicit modeling of survival between embryos and age 0 recruits, and allows the description of a wide range of pre-recruit survival curves. The model is especially useful for low fecundity species that produce relatively few offspring per litter and exhibit a more direct connection between spawning output and recruitment than species generating millions of eggs.

Survival-based recruitment is constrained so that the recruitment rate cannot exceed fecundity. The relationship between survival and spawning output is based on parameters which are on a natural log scale. These are: \[z_0=-ln(S_0)\] which is the negative of the natural log of the equilibrium survival \(S_0\), and can be thought of as a pre-recruit instantaneous mortality rate at equilibrium, and \[z_{\text{min}}=-ln(S_{\text{max}})=z_0(1-z_{\text{frac}})\] which is the negative of the natural log of the maximum pre-recruit survival rate (\(S_{\text{max}}\), the limit as spawning output approaches 0), and is parameterized as a function of \(z_{\text{frac}}\) (which represents the reduction in mortality as a fraction of \(z_0\)) so the expression is well-defined over the parameter range 0 < \(z_{\text{frac}}\) < 1.

Recruitment at age 0 for each year in the time series is calculated as: \[{ R_y = SB_ye^{\Big(-z_0 + (z_0-z_{min})\big(1-(SB_y/SB_0)^\beta \big)\Big)}e^{\tilde{R}_y}\qquad \tilde{R}_y\sim N(0;\sigma^2_R)}\] where \(SB_y\) is the spawning output in year y, \(\beta\) is a parameter controlling the shape of density-dependent relationship between relative spawning depletion \(SB_y/SB_0\) and pre-recruit survival (with limit \(\beta\) < 1), \(\tilde{R}_y\) is the recruitment in year \(y\), and \(\sigma_R\) is the standard deviation of recruitment in natural log space.

As implemented in Stock Synthesis, the parameters needed to apply the stock-recruitment relationship based on the pre-recruit survival are ln(\(R_0\)), \(z_{\text{frac}}\), and \(\beta\). The value of ln(\(R_0\)) defines the equilibrium quantities of \(SB_0\), \(S_0\), and \(z_0\) for a given set of biological inputs (either estimated of fixed), regardless of the values of \(z_{\text{frac}}\) and \(\beta\).

The interpretation of the quantity \(z_0=-ln(S_0)\) as pre-recruit instantaneous mortality rate at unfished equilibrium is imperfect because the recruitment in a given year is calculated as a product of the survival fraction \(S_y\) and the spawning output \(SB_y\) for that same time period so that there is not a 1-year lag between quantification of eggs or pups and recruitment at age 0, which is when recruits are calculated within Stock Synthesis. However, if age 0 or some set of the youngest ages is not selected by any fishery of survey, then density dependent survival may be assumed to occur anywhere before the first appearance of any cohort in the data or model expectations. In such cases, the upper limit on survival up to age \(a\) is given by \(S_{\text{max}}e^{-aM}\).

Nevertheless, interpreting \(z_0\) as an instantaneous mortality helps with the understanding of \(z_{\text{frac}}\). This parameter controls the magnitude of the density-dependent increase in survival associated with a reduction in spawning output. It represents the fraction by which this mortality-like rate is reduced as spawning output is reduced from \(SB_0\) to 0. This is approximately equal to the increase in survival as a fraction of the maximum possible increase in survival. That is: \[z_{\text{frac}}=\frac{ln(S_{\text{max}})-ln(S_0)}{-ln(S_0)} \approx \frac{S_{\text{max}}-S_0}{1-S_0}\] For example, if \(S_0 = 0.4\), \(z_{\text{frac}}=0.8\), then the resulting fraction increase in survival is \((S_{\text{max}}-S_0)/(1-S_0)=0.72\).

The parameter \(\beta\) controls the point where survival changes fastest as a function of spawning depletion. A value of \(\beta\) = 1 corresponds to a linear change in natural log survival and an approximately linear relationship between survival and spawning depletion. Values of \(\beta\)<1 have survival increasing fastest at low spawning output (concave decreasing survival) whereas \(\beta\)>1 has the increase in survival occurring the fastest closer to the unfished equilibrium (convex decreasing survival).

The steepness (\(h\)) of the spawner-recruit curve (defined as recruitment relative to \(R_{0}\) at a spawning depletion level of 0.2) based on pre-recruit survival can be derived from the parameters discussed above according to the relationship and associated inequality: \[h = 0.2e^{z_0z_{\text{frac}}(1-0.2^\beta)}<0.2e^{z_0}=\frac{1}{5S_0}=\frac{SB_0}{5R_0}\]

Unlike the Beverton-Holt stock-recruitment relationship, recruitment can increase above \(R_0\) for stocks that are below \(SB_0\) and thus the steepness is not fundamentally constrained below 1. However, in many cases, steepness will be limited well below 1 by the inequality above, which implies an inverse relationship between the maximum steepness and equilibrium survival. Specifically, the inequality above bounds steepness below 1 for all cases where \(S_0\) > 0.2, which are those with the lowest fecundity, an intuitively reasonable result. For example, with \(S_0\) = 0.4, the steepness is limited below 0.5, regardless of the choice of \(z_{\text{frac}}\) or \(\beta\). This natural limit on steepness may be one of the most valuable aspects of this stock-recruitment relationship.

Code for the survival based recruitment can be found in SS_recruit.tpl, search for “SS_Label_43.3.7 survival based”.

8.6.1.5 Shepherd

The Shepherd stock recruit curve is calculated as: \[R_y = \bigg(\frac{SB_y}{SB_0}\bigg)\frac{5h_{adj}R_0(1-0.2^c)}{(1-5h_{adj}0.2^c)+(5h_{adj}-1)(\frac{SB_y}{SB_0})^c}e^{-0.5b_y\sigma^2_R+\tilde{R}_y}\qquad \tilde{R}_y\sim N(0;\sigma^2_R)\] where \(c\) is the shape parameter for the stock recruitment curve, and \(h_{adj}\) is the transformed steepness parameter defined as: \[h_{adj}=0.2+\bigg(\frac{h-0.2}{0.8}\bigg)\bigg(\frac{1}{5*0.2^c}-0.2\bigg)\]

More details can be found in Punt and Cope (2019).

8.6.1.6 Ricker Re-parameterization

The Ricker stock recruit curve re-parameterized version. More details can be found in Punt and Cope (2019). \[R_y = R_0*(1-temp)*e^{ln(5h)(1-SB_y/SB_0)^{\gamma}/0.8^{\gamma}}\] where \(\gamma\) is the Ricker shape parameter and \(temp\) is defined as: \[temp = \begin{cases} 1-SB_y/SB_0 & \text{if $1-SB_y/SB_0 >$ 0 }\\ 0.001 & \text{if $1-SB_y/SB_0 \leq$ 0} \end{cases}\] where \(temp\) stabilizes recruitment at \(R_0\) if \(SB_y > SB_0\).

8.6.2 Spawner-Recruitment Parameter Setup

Read the required number of long parameter setup lines (e.g., LO, HI, INIT, PRIOR, PRIOR TYPE, SD, PHASE, ..., and BLOCK TYPE). These parameters are:

Example setup of the spawner-recruitment section:

8.6.3 Spawner-Recruitment Time-Varying Parameters

The \(R_{0}\), steepness, and regime shift parameters can be time-varying by blocks, trends, environmental linkages, or random deviations. Details on how to specify time-varying parameters can be found in the Time-Varying Parameter Specification and Setup section. However, not all of these options make sense for all parameters; please see additional details in the section on Time-Varying Stock-Recruitment Considerations.

8.6.4 Tuning \(\sigma_R\)

The \(\sigma_R\) value is typically not estimable and it is recommended practice to tune input \(\sigma_R\) values based on the variance in estimated recruitments post running SS3. The R package R code for Stock Synthesis (r4ss) designed to read and visualize SS3 model results provides recommendations on adjusting \(\sigma_R\) values in the sigma_R_info object in the list created by the r4ss::SS_output() function. An alternative \(\sigma_R\) value is provided based on equation:

\[\sigma_R^2 = Var(\hat{r}) + \overline{SE(\hat{r}_y)}^2\]

Simulation studies conducted by Methot and Taylor (2011) compared three methods of tuning \(\sigma_R\) with the above approach generally performed best since it accounts for variability among recruitment deviations and the uncertainty in their estimates.

8.6.5 Recruitment Deviation Setup

8.6.5.1 Recruitment Eras

Conceptually, the model treats the early, data-poor period, the main data-rich period, and the recent/forecast time period as three eras along a continuum. The user has control of the break year between eras. Each era has its own vector. The early era is defined as a vector (prior to V3.10 this was a deviation vector) so it can have zeros during the earliest years not informed by data and then a few years with non-zero values without imposing a zero-centering on this collection of deviations. The main era can be a vector of simple deviations, or a deviation vector, but it is normally implemented as a deviation vector so that the spawner-recruitment function is its central tendency. The last era does not force a zero-centered deviation vector so that it can have zeros during the actual forecast and non-zero values in last few years of the time series. The early and last eras are optional, but their use can help prevent balancing a preponderance of negative deviations in early years against a preponderance of positive deviations in later years. When the 3 eras are used, it would be typical to turn on the main era during an early model phase, turn on the early era during a later phase, then have the last era turn on in the final phase.

8.6.5.2 Recruitment Likelihood with Bias Adjustment

For each year in the total recruitment deviation time series (early, mid, late/forecast) the contribution of that year to the log likelihood is equal to: \(dev^2/(2.0*\sigma^2_R)+offset*ln(\sigma_R)\) where offset is the recruitment bias adjustment between the arithmetic and geometric mean of expected recruitment for that year. With this approach, years with a zero or small offset value do not contribute to the second component. \(\sigma_R\) may be estimable when there is good data to establish the time series of recruitment deviations.

The implemented recruitment bias adjustment is based upon the work documented in Methot and Taylor (2011) and following the work of Mark N. Maunder and Deriso (2003). The concept is based upon the following logic. The \(\sigma_R\) represents the true variability of recruitment in the population. It provides the constraining penalty for the estimates of recruitment deviations, and it is not affected by data. Where data that are informative about recruitment deviations are available, the total variability in recruitment, \(\sigma_R\), is partitioned into a signal (the variability among the recruitment estimates) and the residual, the variance of each recruitment estimate calculated as:

\[SE(\hat{r}_y)^2 + SD(\hat{r})^2=\Bigg( \bigg( \frac{1}{\sigma^2_d}+\frac{1}{\sigma^2_R}\bigg)^{-1/2}\Bigg)^2+\Bigg( \frac{\sigma^2_R}{(\sigma^2_R+\sigma^2_d)^{1/2}}\Bigg)^2=\sigma^2_R\]

Where there are no data, no signal can be estimated and the individual recruitment deviations collapse towards 0.0 and the variance of each recruitment deviation approaches \(\sigma_R\). Conversely, where there highly informative data about the recruitment deviations, then the variability among the estimated recruitment deviations will approach \(\sigma_R\) and the variance of each recruitment deviation will approach zero. Perfect data will estimate the recruitment time series signal perfectly. Of course, we never have perfect data, so we should always expect the estimated signal (variability among the recruitment deviations) to be less than the true population recruitment variability.

The correct offset (bias adjustment) to apply to the expected value for recruitment is based on the concept that a time series of estimated recruitments should be mean unbiased, not median unbiased, because the biomass of a stock depends upon the cumulative number of recruits, which is dominated by the large recruitments. The degree of offset depends upon the degree of recruitment signal that can be estimated. Where no recruitment signal can be estimated, the median recruitment is the same as the mean recruitment, so no offset is applied. Where log-normal recruitment signal can be estimated, the mean recruitment will be greater than the median recruitment. The value:

\[b_y=\frac{E\Big( SD(\hat{r}_y)\Big)^2}{\sigma^2_R}=1-\frac{SE(\hat{r}_y)^2}{\sigma^2_R}\]

of the offset then depends upon the partitioning of \(\sigma_R\) into between and within recruitment variability. The most appropriate degree of bias adjustment can be approximated from the relationship among \(\sigma_R\), recruitment variability (the signal), and recruitment residual error. Because the quantity and quality of data varies during a time series, the user can control the rate at which the offset is ramped in during the early, data-poor years, and then ramped back to zero for the forecast years. On output, the mean bias adjustment during the early and main eras is calculated, comparing this value to the rmse of estimated recruitment deviations in the Report.sso file. A warning is generated if the rmse is small and the bias adjustment is larger than 2.0 times the ratio of \(RMSE^2\) to \(\sigma^2_R\). Additional information on recruitment bias adjustment can be found in the Recruitment Variability and Bias Correction section.

In mcmc mode, the model still draws recruitment deviations from the log-normal distribution, so the full offset is used such that the expected mean recruitment from this log-normal distribution will stay equal to the mean from the spawner-recruitment curve. When the model reaches the mcmc and MCEVAL phases, all bias adjustment values are set to 1.0 for all active recruitment deviations because the model is now re-sampling from the full log-normal distribution of each recruitment.

8.6.5.3 Recruitment Cycle

When the model is configured such that seasons are modeled as years, the concept of season within year disappears. However, there may be reason to still want to model a repeating pattern in expected recruitment to track an actual seasonal cycle in recruitment. If the recruitment cycle factor is set to a positive integer, this value is interpreted as the number of time units in the cycle and this number of full parameter lines will be read. The cyclic effect is modeled as an exp(p) factor times \(R_{0}, s\)o a parameter value of 0.0 has nil effect. In order to maintain the same number of total recruits over the duration of the cycle, a penalty is introduced so that the cumulative effect of the cycle produces the same number of recruits as \(\text{Ncycles}*R_{0}\). Because the cyclic factor operates as an exponential, this penalty is different from a penalty that would cause the sum of the cyclic factors to be 0.0. This is done by adding a penalty to the parameter likelihood.

8.6.5.4 Initial Age Composition

A non-equilibrium initial age composition is achieved by setting the first year of the recruitment deviations before the model start year. These pre-start year recruitment deviations will be applied to the initial equilibrium age composition to adjust this composition before starting the time series. The model first applies the initial \(F\) level to an equilibrium age composition to get a preliminary N-at-age vector and the catch that comes from applying the \(F\)s to that vector, then it applies the recruitment deviations for the specified number of younger ages in this vector. If the number of estimated ages in the initial age composition is less than maximum age, then the older ages will retain their equilibrium levels. Because the older ages in the initial age composition will have progressively less information from which to estimate their true deviation, the start of the bias adjustment should be set accordingly.

8.7 Fishing Mortality Method

The implementation and reporting of fishing mortality, \(F\), in SS3 is more complex than in simpler models and can be confusing. The description provides an overview of the ways in which \(F\) is calculated, used, and reported.

8.7.0.1 Nomenclature

The nomenclature below ignores sex, morphs and areas for simplicity. The quantities associated with \(F\) calculations are defined as:The quantities associated with \(F\) calculations are defined as:

• \(f\) is fleet.

• \(t\) is a time step; continuous across years \(y\) and seasons \(s\); equivalent to year if only 1 season.

• \(a\) is age.

• \(C_{t,f}\) is fleet-specific catch in a time step.

• \(B_{t,f}\) is fleet specific available biomass, e.g., total biomass filtered by fleet-specific age selectivity, \(s_{t,f,a}\).

• \(s_{t,f,a}\) is age-specific selectivity for a fleet. If selectivity is length-specific, then age-specific selectivity is calculated as the dot product across length bins of length selectivity and the normal (or log-normal) distribution of length-at-age. If selectivity is both length- and age-based, which is an entirely normal concept in SS3, then age selectivity due to length selectivity is calculated first, then multiplied by the direct age selectivity. This compound age selectivity is used in the mortality calculations and is reported as \(Asel2\) in report:32 of Report.sso. See appendix to Methot and Wetzel (2013) for more detail on this.

• \(F_{t,f}'\) is the fully-selected fishing mortality for a fleet. The manual will refer to it as \(F'\) or \(\text{full\_}F\). If your model is using \(F'\text{s}\) as parameters, then the parameter values are for \(F'\).

• \(\text{apical\_}F\) is the maximum \(F\) among the ages. It is not used explicitly in SS3 although some sections of the manual may state \(\text{apical\_}F\) instead of the correct term, full_\(F\).

• \(F_{t,f,a}\) is age and fleet-specific fishing mortality rate equal to \(F_{t,f}' * s_{t,f,a}\). Note that it is possible for no age to have a selectivity equal to 1.0. Some selectivity curves are forced to have maximum selectivity = 1.0 such that \(\text{apical\_}F\) and \(\text{full\_}F\) will be the same, but some curves never reach 1.0 and some non-parametric forms allow selectivity > 1.0. In all cases, \(F'\) is still the rate for the hypothetical age that has selectivity equal to 1.0. The reported \(F'\) values are not rescaled to be an \(F\) for the age with peak selectivity. Users need to take this into account if they are comparing reported \(F'\) values to reported vector of \(F_{t,f,a}\) values.

• \(\text{ann}F_y\) is a measure of the total fishing intensity for a year, based on one of several user-specified options (see below).

• \(F\text{std}_y\) is a standardized measure of the total fishing intensity for a year and is reported in the derived quantities, so variance is calculated for this quantity. See below for how it relates to \(\text{ann}F\).

• \(SPR_y\) is the equilibrium spawning biomass per recruit that would result from fishing at the current year’s \(F\) level and selectivity.

Terminology and reporting of \(\text{ann}F\) and \(F\text{\_std}\) has been slightly revised for clarity in v.3.30.15.00 and the description here follows the new conventions.

8.7.0.2 F_Method

There are two approaches for estimating fishing mortality (\(F\)) in Stock Synthesis. One provides a mid-season harvest rate using Pope’s approximation, the other provides a season-long \(F\) rate using the Baranov catch equation. Pope’s harvest rate is implemented as F_Method = 1. The subsequent options model \(F\) as a season-long rate using the Baranov catch equation.

• F_Method = 2 represents \(F\) as a parameter, \(F_{par}\), for all fishing fleets. Because it is a parameter, its estimation is influenced by all data in the model, with retained catch and discard catch having the greatest influence depending on the standard error for these observations. Starting values for the \(F_{par}\) can be calculated in early model phases using the hybrid approach (below), then copying those \(F_{hyb}\) values into the \(F_{par}\) values beginning in a specified phase for all fleets;

• F_Method = 3 represents \(F\) as a tuned factor, \(F_{hyb}\) using Pope’s and Baranov catch equations sequentially in what is termed the hybrid method. The tuning is based on the retained catch for each fishing fleet;

• F_Method = 4 is a fleet-specific superset of F_Methods 2 and 3 and is the recommended method for modeling \(F\) in Stock Synthesis.

• Discards: SS3 is calculating the \(F'\) to match the retained catch conditional on the fraction of total catch that is retained if the discard/retention feature is used. The \(F\) creates the total catch, then the retention curve partitions it into retained and discarded portions. With hybrid and Pope’s methods, the close tuning to retained catch can result in very imprecise fits to the discard data. If the precision of discard data and catch data are of comparable magnitudes, then it is better to use \(F\) as parameters. SS3 will be using one \(F’\) value to attempt to match two data types. Time-varying retention curves can allow SS3 to fit both catch and discard precisely. But with a non-time-varying retention curve, the fit to catch and discard data will depend on the relative precision of the two data types.

The \(F_{hyb}\) approach works by:

1. Applying the Pope’s method to get a harvest rate for each fleet in the current season;

2. Converting that harvest rate to the equivalent \(F_{hyb}\). This is exact with one fleet and approximate with multiple fleets;

3. Adjusting those \(F_{hyb}\) values over a fixed number of iterations (2-4) using the ratio of observed to calculated catch for each fleet; and

4. Proceeding with those \(F_{hyb}\) values into subsequent model steps. \(F_{par}\) and \(F_{hyb}\) values are used in the same Baranov catch equation and a catch log likelihood is calculated based on the observed retained catch and the annual standard error of the catch.

Some notes:

• At moderate \(F\) levels, hybrid should get a near exact match to the retained catch even if the catch standard error is high. This may result in poor fit to discard catch data if it exists.

• At high \(F\) and/or with many fleets, the hybrid approach may not get to an exact match to the observed catch.

• At high \(F\), the problem becomes quite computationally stiff for Pope’s approximation and the hybrid method so convergence in admb may slow due to more sensitive gradients in the log likelihood with regard to other model parameters. Whereas, admb adjusts the \(F_{par}\) values slowly over the entire sequence of iterations, so it may converge with fewer iterations.

• The hybrid approach performs best when the catch is known with high precision (standard error <   0.05) and the total \(F\) across all fleets is not much greater than natural mortality. The \(F\) as parameter approach works best when total \(F\) is high; also where a fleet has both retained catch and discarded catch, because the \(F\) is influenced by both.

• With F_Method 4, fleets with low \(F\) can apply \(F_{hyb}\) across all model phases and high \(F\) fleets can switch to \(F_{par}\) during later phases.

• It is always advisable to allow the model to start with good starting values for parameters. This can be done by specifying a later phase (> 1) under the conditional input for F_Method = 2 where early phases will use the hybrid method, then switch to \(F\) as parameter in later phases and transfer the hybrid \(F\) values to the parameter initial values. Same advice for using the F_Method = 4. However, bycatch fleets do not have retained catch, so cannot use the hybrid method. They must use \(F_{par}\) with user-specified starting value starting in phase 1.

• Tests have demonstrated that the \(F\) approach has negligible impact on the variance of derived quantities, such as spawning biomass in the final year.

• It is possible to estimate \(F_{par}\) for a particular fleet \(\times\) time even if it has an unknown catch observation (for example, if you have years of catch data followed by a gap and then catch data again, can you estimate the missing catch data in SS3). This certainly is true if there is a discard observation for that fleet \(\times\) time, but the effect of the \(F_{par}\) on the population age composition and trends also allows the survey and composition data to influence the estimates of the \(F_{par}\)s. A catch observation must exist in order to trigger SS3 to estimate a \(F_{par}\) for that fleet \(\times\) time, so enter a dummy (but reasonable) value for the catch observation and give it a very large se so that the deviations from the input catch observations will only be slightly penalized. Make sure that you use F_Method = 4 so that the parameter approach will be used for that fleet. If the \(F\) values for those years fluctuate dramatically, you can stabilize those \(F\) values by adding a dummy effort (survey type = 2) time series. For example, if you are missing catch for 2012–2014, enter an effort value of 1.0 for those years and at least one adjacent year for which SS3 is already able to estimate \(F\). This will stabilize the missing year \(F\)s to be similar to the \(F\)s for the adjacent years.

8.7.1 Equilibrium SPR

\(SPR_y\) is the equilibrium spawning biomass per recruit that would result from fishing at the current year’s \(F\) level and selectivity. It is a measure of the expected long-term effect of fishing. It is calculated as the ratio of the equilibrium reproductive output per recruit that would occur with the current year’s \(F\) intensities and biology, to the equilibrium reproductive output per recruit that would occur with the current year’s biology and no fishing. Thus, it internalizes all seasonality, movement, weird selectivity patterns, and other factors. Because this index moves in the opposite direction than \(F\) itself, it is usually reported as \(1-SPR\). A benefit of this index is that it is a direct measure of common proxies used for \(F_{MSY}\), such as \(F_{40\%}\). A shortcoming of this index is that it does not directly demonstrate the fraction of the stock that is caught each year. The spr value is also calculated in the benchmarks (see below).

The derived quantities report shows an annual spr statistic. The options, as specified in the starter.ss file, are:

• 0 = skip

• 1 = \((1-SPR)/(1-SPR_{TARGET})\)

• 2 = \((1-SPR)/(1-SPR_{MSY})\)

• 3 = \((1-SPR)/(1-SPR_{B_{TARGET}})\)

• 4 = raw spr

8.7.2 Initial Fishing Mortality

Read a short parameter setup line for each fishery and season when there is non-zero equilibrium catch in a season for the fleet (equilibrium catches are input in the catch section of the data file). The parameters are the fishing mortalities for the initial equilibrium catches. Do not try to estimate parameters for fisheries with zero initial equilibrium catch - no parameter line is needed fleets and seasons with zero equilibrium catch.

If there is catch, then give a starting value greater than zero, and it generally is best to estimate the parameter in phase 1. The initial \(F\) parameter lines are ordered as shown in the example below - by season, then within a season, by fleet.

If the initial equilibrium catch is near msy, than a logical inconsistency may occur as documented in the Equilibrium Recruitment section.

It is possible to use the initial \(F\) method to achieve an estimate of the initial equilibrium \(Z\) in cases where the initial equilibrium catch is unknown. To do this requires 2 changes to the input files:

1. Data File: Include a positive value for the initial equilibrium catch for at least one fleet, often with a higher standard error depending upon the situation.

2. Control File: Add an initial \(F\) parameter line (short parameter line) for each fleet and season with initial equilibrium catch to be estimated immediately after the Fishing Mortality setup. It will be influenced by the early age and size comps which should have some information about the early levels of \(Z\).

8.8 Catchability

Catchability is the model process that converts selected numbers or biomass for a fleet into the expected value for a survey or cpue by that fleet. In SS3, the concept has been extended so that, for example, a time series of an environmental factor could be treated as a survey of the time series of deviations for some parameter. Thus, the concept of catchability is expanded to include a family of link functions beyond simple proportionality.

Collectively, all factors that need a catchability are considered to be an index. Note that the catchability options entered in this section of the control file have a strong relationship with the index units (i.e., biomass vs. numbers vs. deviations) and index error type (i.e., normal vs. T-distribution) as entered in the data file.

For each fleet that has index observations, enter a row with six entries as described below. Fleets with no index observations are omitted.

1. Fleet Number

2. Link type or index of deviation vector: An assumed functional form between Q, the expected value, and the survey observation. See example Q parameter lines below for an example set up for link types 4-6.

   1. 1 = simple Q, e.g., proportional: \(y=Q*x\) where \(x\) is a model quantity, like selected biomass, and \(y\) is the expected value for the index.

   2. 2 = mirror simple Q - this will mirror the Q value from another fleet. Mirror in Q must refer to a lower number fleet relative to the fleet with the mirrored Q (example: fleet 3 mirror fleet 2). Requires a Q parameter line for the fleet but will not be used.

   3. 3 = power function to allow non-linearity in the expected value for the index: \(y=Qx^{(1+c)}\). Therefore, \(c > 0\) leads to hyper-stability and \(c < 0\) leads to hyper-depletion. The power function option cannot be used with deviation surveys (type 35 or 36) because of negative values in the deviation vectors.

   4. 4 = mirror Q with scale (2 parameter lines required) with the second parameter (scale) being added to the Q being mirrored. The mirrored Q will be reported as base Q + scale value. Mirror in Q must refer to a lower number fleet relative to the fleet with the mirrored Q. An example usage is when the cpue of a fleet in one area is assumed to have the same catchability as the cpue for a comparable fleet operating in another area. The scale parameter could be used to adjust for the relative spatial areas because Q relates cpue to area-specific population abundance, not population density per area.

   5. 5 = offset (2 parameter lines required). This option, new with v.3.30.23, requires a second parameter that is added to the expected values before multiplying Q. This is most useful for survey types 35 and 36, even if those survey input vectors have already been zero centered. \(y=Q*(x+B)\)

   6. 6 = offset and power (3 parameter lines required). This requires 2 parameters in addition to the Q parameter. The first is for the offset and the second is for the power, \(y=Q*(x+b)^{(1+c)}\).

3. Extra input for link information (i.e., mirror fleet or dev index number)

   1. > 0 = mirror the Q from another (lower numbered survey designated by referencing the fleet number).

   2. If an index of a deviation vector , option 35, use this column to enter the index of the deviation vector to which the index is related. For example, if only using one index in your model, this number would be 1.

   3. If a depletion survey, option 34, approach is being applied the following values in this column determines how phases are adjusted:

      • 0 = add 1 to phases of all parameters. Only \(R_{0}\) active in new phase 1. Mimics the default option of previous model versions;

      • 1 = only \(R_{0}\) active in phase 1. Then finish with no other parameters becoming active; useful for data-limited draws of other fixed parameters. Essentially, this option allows the model to mimic dbsra; and

      • 2 = no phase adjustments, can be used when profiling on fixed \(R_{0}\).

4. Extra_se: Estimable extra standard error for an index

   1. 0 = skip (typical); and

   2. 1 = include a parameter that will contain a value to be added to the input standard deviation of the survey variability.

5. Bias_adj: Adjusts for log-normal bias when using an informative prior on Q. NOTE: Bias adjustment to Q is ONLY applied when also using the float approach. See below for Q float specifications. An expanded bias adjustment approach is under development.

   1. 0 = no bias adjustment applied; and

   2. 1 = apply bias adjustment. Bias correction will be applied to the estimated Q value.

6. Float: Allows Q parameter to be automatically calculated to maintain no deviation between the index and the expected value. Only available for link = 1.

   1. 0 = no float, parameter is estimated; and

   2. 1 = float, analytical solution is used, but parameter line still required.

   3. Additional information regarding the use of and appropriate application of float in Q can be found in the Float Q section below.

If the Q base parameter specifies that it is time-varying by the annual deviation method, short parameter lines to specify the specifications of the deviation vector come after all the base Q parameters.

8.8.1 Example Q parameter lines

Below is an example setup for fleets with a mirrored Q and scale from another fleet (link type = 4), scale from another fleet with offset (link type = 5), and scale from another fleet with offset and power (link type = 5):

8.8.2 Float Q

The use and development of float in Q has evolved over time within SS3. The original approach in earlier versions of SS3 (version 3.24 and before) is that with Q “float” the units of the survey or fishery cpue were treated as dimensionless so the Q was adjusted within each model iteration to maintain a mean difference of 0.0 between observed and expected (usually in natural log space). In contrast, Q as a parameter (float = 0) one had the ability to interpret the absolute scaling of Q and put a prior on it to help guide the model solution. Also, with Q as a parameter the code allowed for Q to be time-varying.

Then midway through the evolution of the v.3.24 code lineage a new Q option was introduced based on user recommendations. This option allowed Q to float and to compare the resulting Q value to a prior, hence the information in that prior would pull the model solution in direction of a floated Q that came close to the prior.

Currently, in v.3.30, that float with prior capability is fully embraced. All fleets that have any survey or cpue options need to have a catchability specification and get a base Q parameter in the list. Any of these Q’s can be either:

• Fixed: by not floating and not estimating.

• Floating: not estimated as an active parameter, set phase to negative, and not capable of being time-varying. Can have a prior, or not. Future versions may allow Q to be time-varying and then rescaling that Q vector according to the float logic, but this is not yet currently implemented.

• Estimated as active parameter: so not floating. Can be time-varying and can have a prior.

Q relates the units of the survey or cpue to the population abundance, not the population density per unit area. But many surveys and most fishery cpue is a proportional to mean fish density per unit area. This does not have any impact in a one area model because the role of area is absorbed into the value of Q. In a multi-area model, one may want to assert that the relative difference in cpue between two areas is informative about the relative abundance between the areas. However, cpue is a measure of fish density per unit area, so one may want to multiply cpue by area before putting the data into the model so that asserting the same Q for the two areas will be informative about relative abundance.

In v.3.30.13, a new catchability option has been added that allows Q to be mirrored and to add an offset to ln(Q) of the primary area when calculating the ln(Q) for the dependent area. The offset is a parameter and, hence, can be estimated and have a prior. This option allows the cpue data to stay in density units and the effect of relative stock area is contained in the value of the ln(Q) offset.

8.8.3 Catchability Time-Varying Parameters

Time-Varying catchability can be used. Details on how to specify time-varying parameters can be found in the Time-Varying Parameter Specification and Setup section.

8.8.4 Q Conversion Issues Between Stock Synthesis v.3.24 and v.3.30

In v.3.24 it was common to use the deviation approach implemented as if it was survey specific blocks to create a time-varying Q for a single survey. In some cases, only one year’s deviation was made active in order to implement, in effect, a block for Q. The transition executable (ss_trans.exe) cannot convert this, but an analogous approach is available in v.3.30 because true blocks can now be used, as well as environmental links and annual deviations. Also note that deviations in v.3.24 were survey specific (so no parameter for years with no survey). In v.3.30, deviations are always year-specific, so you might have a deviation created for a year with no survey.

8.9 Selectivity and Discard

For each fleet and survey, read a definition line for size selectivity and retention.