Sensitivity Analysis in a Complex Marine Ecological Model

Sensitivity analysis (SA) has long been recognized as part of best practices to assess if any particular model can be suitable to inform decisions, despite its uncertainties. SA is a commonly used approach for identifying important parameters that dominate model behavior. As such, SA address two elementary questions in the modeling exercise, namely, how sensitive is the model to changes in individual parameter values, and which parameters or associated processes have more influence on the results. In this paper we report on a local SA performed on a complex marine biogeochemical model that simulates oxygen, organic matter and nutrient cycles (N, P and Si) in the water column, and well as the dynamics of biological groups such as producers, consumers and decomposers. SA was performed using a “one at a time” parameter perturbation method, and a color-code matrix was developed for result visualization. The outcome of this study was the identification of key parameters influencing model performance, a particularly helpful insight for the subsequent calibration exercise. Also, the color-code matrix methodology proved to be effective for a clear identification of the parameters with most impact on selected variables of the model.


Introduction
Sensitivity analysis (SA) can be basically described as the process to evaluate the contribution of input parameters to model behavior.Several parameters in ecosystem models represent specific process coefficients that are only measured with difficulty (if possible at all).As a consequence, there are uncertainties related with the parameterization and the nonlinearity of interactions within the model.This raises two basic questions: (a) how sensitive is the model to changes in individual parameter values; and (b) which parameters or associated processes have more influence on specific output variables?Answering the first question may reveal which parameters have more influence on the model and need more attention.The answer to the second question will help to understand the simulated system.
Performing a SA is an effective way to answer these questions and assess model performance [1][2][3][4][5], thus justifying its widespread use to quantitatively assess the influence of parameters in model performance, and identify those having the most significant impact on the results [6][7][8][9][10][11][12].In complex models with a significant number of parameters, the SA may help in the selection of the most relevant parameters for the calibration process [4,7,13,14].
Sometimes the large number of input values in the habitual "one at a time" SA perturbation method requires excessive computation time [14][15][16].Also, it leads to serious difficulties in visualizing the results in a comprehensive manner, especially in complex ecosystem models with dozens of parameters.To avoid a sensitivity matrix with numerous columns or rows, some proposed methodologies do not analyze single parameters, but rather clusters of related parameters in group-collecting sensitivity analysis, reducing the total amount of required variations by at least a factor 1/9 [17].Other approaches, such as metamodeling, link model inputs to the output through a known relationship, emulating the original model but with less computational demands [10,18].Despite all the advantages of such SA tests, one of the drawbacks is that it can mask the importance of individual parameters in model performance, providing instead a black-box response for the effect of individual parameters [19].
The aim of this paper is to perform a SA on a complex marine biogeochemical model that reproduces the dynamics of living groups (producers, consumers and decomposers), main nutrients (N, P and Si), oxygen and several pools of organic matter.We have used a schematic mesocosm application to prevent any influence of physical transport processes on the outcome of the model.To simplify the analysis, we have developed a color-code matrix to visualize the results of the SA.With this study we intend to identify key parameters influencing model performance, to focus on their variation in the subsequent calibration process.

Model and Methods
A process-oriented zero-dimensional ecological model (MOHID-Life [20]) is applied to simulate biogeochemical constituents in a virtual mesocosm.The choice of a zero-dimensional model application, instead of a real case application, is justified because it avoids a substantial increase in model run time.Also, it guarantees that only ecological processes are responsible for the change in the state-variable over time, by eliminating the effect of transport on the state-variables.The only physical process considered in the simulations is the sinking of particulate organic matter.MOHID-Life consists of a series of coupled first-order differential equations representing the major biogeochemical processes influencing the water quality, as well as the dynamic of several groups of primary and secondary producers, the microbial loop, nutrient recycling and oxygen dynamics.MOHID-Life is a complex model that accounts for the cycles of carbon, nitrogen, phosphorus, silica and oxygen.A detailed description of the model options and philosophy, as well as its parameterization, has been previously described and the values of the parameters tested here have been derived from these studies [20,21].
The flexible parameterization of the model allows configuration of different degrees of complexity, by reducing or increasing the number of phytoplankton and zooplankton groups.So, to reduce the volume of results generated by the runs, and aiming at simplicity in the analysis, the reference simulation for this study is simplified by considering only one producer (diatoms) and one consumer (microzooplankton).Because the model philosophy is built on the concept of a Generic Type Model [20], this simplification in the application still allows assessing the importance of parameters, since the model uses the same code for all groups of producers and for consumer.
A reference run is used to provide the benchmark results for the parameter sensitivity analysis (see Section 2.3).The results of the sensitivity analysis are classified by distinguishing model parameters with a qualitatively different effect on the outcome of the simulation.The effect of single parameter variation on model performance is analyzed by singly and sequentially altering the standard value of each parameter with up-and down-variations of 10% in a series of separate runs, while holding all other terms constant.A raised and lowered 10% parameter perturbation is frequently adopted [15,22,23], hence its use in the present study.

Sensitivity Index
The choice of a sensitive index varies significantly in the literature [1,10,24,25].From all the indexes available to quantify parameter sensitivity, we chose one able to quantify normalized sensitivity, but with the capacity to reveal up (positive) and down (negative) variation.The choice of index was based on its simplicity, and directly related to the objective of knowing the influence of any given parameter on the outcome of the simulation, i.e., the impact of a parameter and the nature of variation in the results.
Normalized sensitivity, S(p), is defined as the relative change in model output divided by the relative change in the parameter value.It is calculated as: where all variables with an S in the lower index represent standard case values (Vs the value of a given variable for the standard case with parameter ps), and V(p) is the value for the case when the parameter is given the value p.As an example, to evaluate the impact of a negative perturbation (−10%) of the parameter bio_si_diss on silicate concentration, we would have ps = 0.01 and p = 0.009, and the respective model results for each parameter value, Vs = 3.07 and V(p) = 2.85, respectively.In this particular case, S(p) = 0.071.This method is proposed by Fasham et al. [26] and adopted in other studies (e.g., [23,27]).According to this sensitivity index, a negative parameter perturbation (10% below) with a negative index result means a positive perturbation in the end result of a given property (meaning a higher end value compared with the reference run value), whereas a positive index result means a negative perturbation in the result (a lower end value compared with the reference run).Conversely, a positive parameter perturbation (10% above) will give a negative index result if a negative end result is achieved, and a positive index result with a positive end result.
The degree of model sensitivity towards any given parameters can be defined as sensitive (S > 0.1, meaning a change of more than 1% in the result when compared with the reference value), highly sensitive (S > 1, meaning a change of more than 10%), and extremely sensitive (S > 10, meaning a change of more than 100%).Whenever S < 0.1, it can be said that the model is not sensitive to that parameter.However, it must be kept in mind that the bias achieved by omitting some variations is a frequent or potential error that might occur in a systematic SA.In addition, it is difficult to define the magnitude of parameter perturbation avoiding non-realistic values, but at the same time covering the actual range of the parameter [26].

Variables of Interest
Not all state-variables were selected to monitor the sensitivity of parameters.This choice relied on the assumption that an excessive number of state-variables would compromise the comprehensibility of the analysis.The biomass of functional groups, concentration of nutrients and labile-DOC (DOMl), and chlorophyll content in phytoplankton were selected as sensitivity indicators (Table 1), mostly because they correspond to properties for which field data is more frequently available.The phytoplankton carbon content to bacterioplankton carbon content ratio (Pc:Bc) was set to evaluate the varying microbial community composition in response to different parameters values.
Table 1.Selected variables to assess model sensitivity.

Symbol Description Initial Value P c
Producers biomass Producers-decomposer ratio 1

Simulations Runs
The model was applied to a virtual mesocosm corresponding to a schematic reservoir with 3 × 3 × 5 square cells with 1 m each.A spin-up period of 16 months was used to stabilize the model, after which the SA simulations run for 3 months (May to July).The long spin-up period aimed at stabilizing the model after the "chaotic" oscillations characteristic from the first year run.The standard model parameterization used in the simulations is presented in Table 2 to Table 5, and have been taken from previous studies [20,21].The initial values for the properties addressed in the SA are provided in Table 1.Environmental conditions used in the model forcing are shown in Figure 1.A schematic water temperature time-series was used and solar radiation corresponds to conditions representative for mid-latitudes in the northern hemisphere.The term "variable" is loosely used throughout the sensitivity analysis discussion because, with the exception of Pc:Bc, the remaining correspond to state-variables in model.

Parameter Description Reference Value Units
Affinity for NH 4 (uptake rate) 0.0025 Affinity for NO 3 (uptake rate) 0.0025 Affinity for PO 4 (uptake rate) 0.0025  Affinity for PO 4 (uptake rate) 0.0025 Affinity for NH 4 (uptake rate) 0.0025 Affinity for NO 3 (uptake rate) 0.0025

Interpretation of the Results
A simple Visual Basic (VBA) macro was developed in Excel environment to process the huge amount of data produced by the model outputs for all the simulations.The macro loads the model results into a spreadsheet and then calculates the sensitivity index of the selected variables to each parameter.In the next step the matrix is interpreted to convert the quantitative index into a qualitative descriptor (the degree of model sensitivity), attributing the + and − symbols to indicate a positive or negative perturbation, respectively.The degree of sensitivity is then given by + and − for sensitive (S > 0.1), ++ or −− for highly sensitive (S > 1), and + + + or − − − extremely sensitive (S > 10).Finally, to allow a better visualization of the results, a color code is used to identify the deferent degrees of sensitivity, positive or negative.No symbol or color is assigned when the model is not sensitive to a given parameter (S < 0.1).A step-wise description of this process is presented in Tables 6-8.Table 6.Representation of the processing and interpretation of the SA result, Step 1: the sensitivity index is calculated from the results of the simulations.Only a fraction of the results are shown for illustrative purposes.The description of parameters and variables can be found in Tables 1-5.

Results
The effect of 68 parameters was tested for their positive and negative variation influence on 11 variables, resulting in a matrix of results with 1496 values.For simplicity, only the most relevant parameter influences on model results, i.e., denoting moderate or high model sensitivity toward that parameter, are discussed.The results are presented in tables for each functional group parameters and for the general parameters, where the impact of each parameter perturbation is expressed in a qualitative way.The tables provide a condensed view of the complexity of the interrelations between parameters and variables within the model.
The simulated period was chosen to have temporal variability, namely the formation and destruction of spring bloom of phytoplankton, thus avoiding a period of model stasis.Results of the reference run are illustrated in Figure 2.

General Parameters
Variables sensitivity to general parameters is represented in Table 9, where it is possible to see that only three parameters (bio_si_diss, pom_bac_ks and pom_bac_vmax) have a significant impact on the results (>1%).Biogenic silica dissolution rate (bio_si_diss) perturbation induces a moderate change in results, especially in producers (also reflected in chlorophyll) and consequently on the Pc:Bc ratio.
Pom_bac_ks and pom_bac_vmax induce small changes in nutrient variables, but have a high impact on producers (and all related variables), especially pom_bac_vmax (S > 10%).A negative perturbation of this parameter leads to the increase in the biomass of producers, and a positive perturbation has the contrary effect.This parameter controls the rate at which hydrolysis converts POM to DOM, a substrate for decomposers.This means the higher pom_bac_vmax is, more DOM is made available for decomposers; these will compete with producers for nutrients thus impacting their growth.

Producers' Parameters
The SA of producers' parameters (Table 10), shows a minor impact of their perturbation in the result of consumers and decomposers, the exception being for max_assimil and ref_temp with a sensitivity of S > 0.1 in consumers.Overall, nutrient affinity to both forms of nitrogen were the parameters with less effect on any of the variables.The model is extremely sensitive (S > 10) only to three parameters, namely, exudation under nutrient stress (exu_nut_stress), maximum assimilation rate (max_assimil), and reference temperature (ref_temp).The effect of these parameters is particularly strong in all phytoplankton related variables (biomass and Chla content), a fact also observed in other studies [15].In this particular case, the simulation ends in a situation of apparent nutrient shortage, and as nutrient limitation increases, the influence of parameters governing growth and response to nutrient stress increases.This is observed in the decrease of exu_nut_stress, meaning less exudation under nutrient stress, resulting in a positive perturbation on the biomass of producers.
The highest observed sensitivity is for the ref_temp upper perturbation in Chla.Overall, ref_temp is the parameter to which producers are more sensitive, an expected occurrence given the control of temperature on several physiological processes.

Consumers' Parameters
Globally, the model is more sensitive to consumers' parameters.SA results in Table 11 show that only some variables are sensitive (S > 0.1) to minimum and maximum nutrient:carbon ratios, and that the model is highly or extremely sensitive to all other parameter, especially the variables related to phytoplankton (biomass, Chla, Pc:Bc).Consumers parameterization is important to the control of both producers and consumers groups, and not so relevant for the decomposers and organic matter dynamics.
The highest sensitivity values observed are for assimil_effic and p_graz_avail by Chla and Pc.The sensitive index observed in Chla for assimil_effic has the highest value in the SA. Figure 3 illustrates the influence of these parameters on the temporal evolution of some variables, and it is possible to notice that, while Zc and Bc are not as sensitive as Pc, their dynamics is affected by this parameter.When compared to the parameter sets of other groups, consumer parameters have the strongest effect on model behavior, most probably due to the grazing control that consumers have on producers.Table 9. SA result matrix for the general parameters impact on the analyzed variables.Color code: blue for sensitive (S > 0.1), orange for highly sensitive (S > 1) and red for extremely sensitive (S > 10).The signs + and − stand for positive and negative perturbation, respectively.A description of the variables and parameters is found in Tables 1-5. Parameter Table 10.SA result matrix for producers' parameters impact on the analyzed variables.Color code: blue for sensitive (S > 0.1), orange for highly sensitive (S > 1) and red for extremely sensitive (S > 10).The signs + and − stand for positive and negative perturbation, respectively.A description of the variables and parameters is found in Tables 1 and 2 11.SA result matrix for consumers' parameters impact on the analyzed variables.Color code: blue for sensitive (S > 0.1), orange for highly sensitive (S > 1) and red for extremely sensitive (S > 10).The signs + and − stand for positive and negative perturbation, respectively.A description of the variables and parameters is found in Tables 1 and 3, respectively.
Table 12.SA result matrix for decomposers parameters impact on the analyzed variables.Color code: blue for sensitive (S > 0.1), orange for highly sensitive (S > 1) and red for extremely sensitive (S > 10).The signs + and − stand for positive and negative perturbation, respectively.A description of the variables and parameters is found in Tables 1 and 4, respectively.

Decomposers' Parameters
The results for the SA analysis on decomposers parameters are illustrated in Table 12.The model only shows to be extremely sensitive (S > 10) to the positive perturbation of ref_temp, and the effect of the perturbation in the outcome of the model can be seen in Figure 4. Despite the oscillation seen in the time series, the simulations converge to a similar result.As in producers, the strong control of temperature on physiologic functions makes the model extremely sensitivity to temperature reference values.Producers and related variables are the most sensitive to decomposers parameter perturbation.This is an expected outcome, since the microbial loop control much of the available nutrients (via respiration of organic matter), affecting directly the dynamics of nutrient uptake by producers.This effect, however, is not noticed in the nutrient variables, meaning that they are being consumed as they become available.DOMl also shows sensitivity to the perturbation of the parameters that control its mineralization by bacteria, as seen in Figure 4.The affinity for nutrient does not show any relevant disturbance on model response, in a similar way as observed for producers.

Discussion
Generally, the model seems to be sensitive to most parameters.The observed sensitivity, however, is not detected in all model compartments, implying that a single parameter may have an impact on one or several state-variables, but not on the overall outcome of the model.Similar observations have been reported in similar studies [3,4,25,28], and can be explained by the complexity of the model in which all state-variables are dependent on processes controlled by multiple parameters.
The model is extremely sensitive to only three parameters in producers: Exudation under nutrient stress (nut_stress_thresh), maximum assimilation rate (max_assimil), and reference temperature.Overall, reference temperature (ref_temp) is the parameter to which producers are more sensitive, an expected occurrence given the control of temperature on several physiological processes [28,29].Results show that consumer's parameterization is particularly relevant to the control of producers, and not so important for decomposers.This is observed in the number of parameters in consumers for which Pc and Chla show high or extremely high sensitivity (Table 11).These state-variables are particularly sensitive to parameters related with the feeding and assimilation in consumers (assimil_effic, max_spec_up_@10c, graz_up_ks and p_graz_avail).
Consumers' parameters have the strongest effect on the behavior of the model, probably through grazing control on producers.Some studies have shown a similar control of consumers (zooplankton) parameters on model performance [2,15,30,31], denoting the importance of this group in ecological models for pelagic systems.Decomposers have the lowest number of parameters for which variables are extremely sensitive.This is only observed for the reference temperature (ref_temp) for this functional group (Figure 4).
Of all variables studied in this SA, nutrients were the least affected in the number of parameters to which they are sensitive.Nutrients sensitivity to parameter perturbation never exceeded 10% in the final result.From this, it is possible to conclude that nutrient variables are the most robust.It can be hypothesized, however, that model sensitivity to parameters related with nutrient dynamics may increase under nutrient limitation.
Silica is the variable affected by the lower number of different parameters, because it is the only variable depending only on producers.Remarkably, no producer parameter disruption causes any significant change in the final outcome of silica.Biogenic silica dissolution rate (bio_si_diss) perturbation produced a minor impact (<10%) on silica, implying that the model is not very sensitive to it.Nevertheless, this parameter showed to be highly sensitive to other state-variables, namely Pc and Chla.This impact is probably due to limitation by silica in producers.Small variations in these parameters will affect the amount of dissolved silica to diatoms, because it mediates the conversion of biogenic silica to dissolved silica.Also, the other parameters to which the model silica compartment is sensitive are related to consumers grazing activity (ref_temp, assimil_effic and max_spec_up_@10c).The grazing pressure influences silica abundance, because this element is not assimilated by consumers, but instead diverted to the biogenic silica pool.
Phosphate only reflects the variation of two parameters, the upper perturbation of max_pc_ratio and the lower perturbation of min_nc_ratio (Table 11), because they disturb the immobilization and mineralization of phosphorus.Hence, their perturbation has a direct impact on phosphate, even though not a significant impact (<10%).Nitrate is not a product of mineralization, so it is not affected by any consumers or decomposer parameters, and is only affected by the nitrification rate perturbation (both positive and negative) and by max_nc_ratio.Of all groups, consumers have more parameters affecting ammonium dynamics, probably because of its important role in nutrient mineralization through grazing and subsequent nutrient recycling.Most of the parameters related to grazing activity have an effect on ammonium.
Ultimately, one of the main outcomes of SA is to identify where a relatively minor change in a factor's value leads to a major change in model output [32,33].This is an essential procedure in conventional model development, illustrated by the volume of published works on the topic and stressed in official guidelines regarding the modeling of environmental systems (see [34,35]).In our study we opted for a local SA, where only one parameter varies at a time, while all the others remain fixed at predetermined values, instead of a global analysis in which all parameters vary [7,36].In global SA, or total SA (sensu [37]), all factors are changed together across the full multidimensional input space.This means that the main effect of parameters and the interactions terms involving them are both included, enabling to determine if parameter influence is mostly due to interactions terms rather than to linear effects.Several studies on complex ecosystem models have shown that global SA are more useful because they address parameter interactions and potential nonlinear relationships between parameters and model outputs [12,[38][39][40][41].
With the simple parameter perturbation approach used in this study it is only possible to suppose potential interactions between parameters controlling model behavior.Taking the most sensitive parameter in consumers as an example, ref_temp and max_spec_up_@10c, it is possible to infer a potential interaction in the control of consumers, producers and nutrients, because both parameters are used to determine the grazing rate.The same could be inferred for two of the most sensitive parameters in producers, ref_temp and max_assimil, since both mediate the maximum rate of photosynthesis [20].
Global SA may provide useful information on the most uncertain parameters, and it certainly shows that model parameterization should be done with carefully determined empirical values because of the interaction effects [29].Nonetheless, a drawback in this approach relies on the extensive computation needed to address all interaction terms and total effects for each input parameter [41,42].Because global SA aims to characterize models performance and uncertainty [42], it may represent an excessive effort when the aim is to identify key parameters to address in the calibration of the model.Similarly, it may not be as effective as local SA on illustrating the direct relationship between the perturbations of a single parameter and model output, given the volume of results produced.While not providing information on nonlinear responses and parameter interaction (or their combined effect), local SA still provides fundamental information to identify the most sensitive parameters governing model behavior.This is significant information because the factors the model is most sensitive to usually have stronger main rather than interaction effects [29].
Graphical representation of SA results, especially global SA, has become a challenge of its own given the continuously increase in complexity of ecological models.A way to reduce the volume of results, and to achieve simplification, relies on grouping parameters and assessing model sensitivity to these assemblages [12,[43][44][45].While these representations have its own benefits, they also have the obvious disadvantage of covering-up the impact of single parameters.This is, in fact, one of the aims of this methodology [46][47][48][49].Since the objective of the work presented here was to identify the key parameters controlling the model behavior, it relied on an approach that explicitly addresses the impact of single parameters.This is a standard and essential practice in process-oriented water-quality and ecological models, because it precedes the calibration exercise [4,5,[50][51][52][53][54][55].Despite the argumentation on the constraints of the one-factor-at-a-time SA [4,6,7,28], we have simplified the complexity in the interpretation of results by using a color-code matrix, thus minimizing some of its drawbacks.With this approaches it becomes easier to relate parameters (model input) and state-variables (model output) and have instantaneous perception on the magnitude of this relationship.When limited data sets are available (for a few state-variables, for instance), the color-code tables may well be an efficient way to calibrate a model adjusting only a few parameters.

Conclusions
Most parameters exert some influence on the selected outputs of the model, but the study shows that the proposed degree of variation in the standard parameterization does not yield a much different scenario, or leads to any significant change in the model results.In the proposed methodology, the result tables provide a condensed view on the sensitivity of state-variables to individual parameters.These tables simplify the interpretation of SA methods that produce significant amounts of results, and can be used as look-up tables to identify sensitive responses in real cases applications.
The methodology presented here is based on local SA.While not providing the same volume of information on the model behavior as a global SA, it yields fundamental information to identify the most sensitive parameters governing model behavior.As such, our approach proved to be a simple and objective way to assist in the calibration exercise of complex models.

Figure 1 .
Figure 1.Environmental conditions used to force the model: (a) Surface temperature and (b) Daily maximum solar radiation.

Figure 3 .
Figure 3. Temporal evolution of the response of (a) producers biomass; (b) chlorophyll a (Chla); (c) consumers biomass; and (d) decomposers biomass to the perturbation of assimil_effic.

Table 5 .
General parameters (not necessarily related with any functional group).

Table 7 .
Representation of the processing and interpretation of the SA result, Step 2: the sensitivity index is arranged in a matrix that relates input parameters and output variables.

Table 8 .
Representation of the processing and interpretation of the SA result, Step 3: for a straightforward interpretation, the sensitivity value is converted to a qualitative descriptor, using signs (+ and −) and a color scale to indicate positive or negative perturbation, and the degree of sensitivity of a variable to a given parameter.