How morphology shapes the parameter sensitivity of lake ecosystem models

A global sensitivity analysis of a lake ecosystem model (GOTM-FABM-PCLake) was undertaken to test the impacts of lake morphology on parameter sensitivity in three different lakes. The analysis was facilitated by the Parallel Sensitivity and Auto-Calibration tool (parsac) and included a screening step with the density-based Borgonovo ’ s method followed by in-depth analysis with both Borgonovo ’ s and the variance-based Sobol ’ methods. The Borgonovo ’ s method proved efficient in ranking the most influential parameters and its results were corroborated by the Sobol ’ method. For total phosphorus and total nitrogen, parameters related to the benthic-pelagic coupling and phytoplankton, were particularly important for the shallower lakes, whereas the most important parameters for total nitrogen were related mainly to the benthic-pelagic coupling in the deepest lake. For chlorophyll a , phytoplankton and zooplankton parameters were most influential. We conclude that lake morphology shapes the parameter sensitivity of lake ecosystem models.


Introduction
Lake ecosystems and their biodiversity are currently under serious threat from eutrophication and global climate changes (Albert et al., 2020;Jeppesen et al., 2014;WWF, 2018).Therefore, lake ecosystem models are increasingly applied to simulate water quality and to gain a deeper understanding of lake ecosystem functioning (Trolle et al., 2012).Lake ecosystem models may also act as tools supporting management decisions and policy making, for instance by quantifying the effectiveness of management measures (Hilt et al., 2017), thereby providing information on how to best mitigate future impacts (e.g.Rolighed et al., 2016).Several complex aquatic ecosystem models are available (Mooij et al., 2010) and have been used to assess the effect of different management measures and the impacts of future global climate changes.For example, (i) three different process-based lake ecosystem models were used to project the response of phytoplankton to climate change in a small, shallow Danish lake (Trolle et al., 2014); (ii) ensemble modelling with a hydrodynamic lake model explored the impacts of climate warming and increased frequency of heat waves on the thermal characteristics of a large sub-tropical lake (Gal et al., 2020) and (iii) the lake ecosystem model complex GOTM-FABM-PCLake was used to assess the level of external nutrient load reductions needed for a shallow lake to reach clear water conditions with a high coverage of submerged vegetation (Andersen et al., 2020).
The algorithms of complex environmental models often include hundreds of parameters, many of which require tailored configuration and adjustment in each application case, and model complexity will likely increase further in the future (Wagener and Pianosi, 2019).This is also true for lake ecosystem models (Robson, 2014), and the majority of model parameters cannot be accurately measured or are unknown due to, for instance, lack of observations.This means that modellers have to rely on literature values and calibration to identify the optimal parameter values to give a satisfactory description of the system being modelled.To ensure robust and high quality model output, it is vitally important to focus on accurately determining the parameters that are most influential for the model outputs (Arhonditsis and Brett, 2004).Comprehensive sensitivity analysis (SA) can serve this purpose and is considered an important part of any due diligent modelling process (Jørgensen, 1995;Ravalico et al., 2005), especially in a policy context (Saltelli and Funtowicz, 2014).This is also stressed by several agencies that recommend the use of SA as best practice in the development and application of models (European Commission, 2005; U.S. EPA, 2009).
SA methods have commonly been divided into two broad categories: local SA methods, where input parameters has been varied one at a time around a nominal value and global SA methods which attempts to explore the entire space of input parameters.Global SA is to be preferred in most cases as (i) it does not assume that the model is linear or additive, which is rarely the case for any complex ecological model, and because (ii) it accounts for the effects of interactions between model parameters (Saltelli et al., 2008;Saltelli and Annoni, 2010).Local SA is still the most commonly used approach, but global approaches have recently achieved increasing attention (Ferretti et al., 2016).
In the development and application of the most widely used lake models, SA has usually been performed and reported, but SA is typically not an integral and ongoing process in further model developments (Mooij et al., 2010).Application of lake models to specific case studies tends to rely on SA results to decide which parameters and their ranges to include in the calibration process (e.g.Di Maggio et al., 2016;Rolighed et al., 2016).SA of lake models has either been undertaken using local (e.g.Schladow and Hamilton, 1997;Elliott et al., 1999;Nielsen et al., 2014) or global (e.g. Bruce et al., 2006;Janse et al., 2010;Page et al., 2017) SA methods.Due to the complexity of ecological models, SA computation can be time-consuming and computational costly, and time and resources are only rarely available for conducting rigorous global SA of such models (Jakeman et al., 2006;Robson et al., 2008).
The sensitivity and calibration utility parsac (Parallel Sensitivity Analysis and Calibration) has the potential to greatly reduce the time required for conducting global SA by automating key steps in the procedure.parsac is written in Python for models with output in NetCDF and model configuration files in namelist or YAML format.To further reduce the computational burden, the global SA may be divided into a screening step and a more detailed analysis step, as recommended by Saltelli et al. (2008).The screening step aims to identify all non-influential parameters conditional on the chosen target variable and concentrate the detailed global SA analysis on influential parameters.This approach has also been applied to lake models with various SA methods (e.g.Janse et al., 2010;Makler-Pick et al., 2011).
Parameterisation of environmental models can differ between sites and contexts, and parameter sensitivity has also been shown to be potentially context specific (Wagener and Pianosi, 2019).For example, parameter sensitivity of the Soil and Water Assessment Tool (SWAT) varied in two different watersheds and with different climatic conditions (Cibin et al., 2010).This is likely also true for lake ecosystem models.Even though most lake models have been developed to cover the general behaviour and mechanisms of lake ecosystems, they often need to be calibrated for a specific lake.Consequently, the analysis may only provide information on influential parameters and model behaviour for that specific lake or lake type.
In this study, we aim to evaluate the influence of lake morphology on the sensitivity of model parameters.We provide the first comprehensive global sensitivity study of the 1-dimensional hydrodynamic-lake ecosystem model complex GOTM-FABM-PCLake for three Danish lakes with varying bathymetry (max.depth ranging from 2 to 33 m) and trophic state (ranging from meso-to eutrophic).Model applications for these three lakes have previously undergone rigorous calibration and validation (Andersen et al., 2020;Chen et al., 2019;Trolle et al., 2008b).As shallow lakes likely have a tighter coupling between the pelagic-benthic domains compared with deep lakes (Schindler and Scheuerell, 2002), we expect to identify more influential model parameters relating to this particular coupling in the shallow lakes.
SA has proven to be sensitive to the choice of methods (Wagener and Pianosi, 2019).We therefore applied two global SA methods: the density-based Borgonovo's δ-sensitivity measure and the widely-used variance-based first-and total-order Sobol' indices in a two-step approach with a screening followed by an in-depth analysis.The global SA was performed with the sensitivity and calibration utility parsac version 0.5.7, which was chosen for its ability to minimise the computational burden and included both Borgonovo's and the Sobol' methods.
The aim of the sensitivity analysis was to obtain insight into which model parameters are most influential in shaping model behaviour across different lake depths.We provide a ranking of the most influential parameters for each lake on important model outputs, which are typically targeted for calibration purposes.We also exemplify usage of the sensitivity and calibration utility tool parsac.The insight from sensitivity analysis in this work will be broadly applicable in future lake model case studies within the temperate zone and for the development of the next generation of lake ecosystem models.

The study design
The focus of this sensitivity analysis is on the output variables that are most commonly examined in calibration and validation processes during a lake modelling case study (Table 1).The analysis includes three lakes with varying bathymetry and trophic state, all located in Denmark: Lake Ravn, Lake Bryrup and Lake Hinge.For each lake, the GOTM-FABM-PCLake model complex have been set up, calibrated and validated in separate studies (for details see Andersen et al. (2020); Chen et al. (2019);Trolle et al. (2008a)).To minimise effects of model initialisation, model simulations included a warm-up period of 5 years (1996)(1997)(1998)(1999)(2000) before the analysis period of 5 years (2001)(2002)(2003)(2004)(2005) and were configured with a 1-hourly time step and daily output values for layers in the surface (ranging between 0 and 3 m).All forcing and boundary conditions were lake specific.Output variables were processed across the analysis period into summer (April to September) mean values.
All 356 FABM-PCLake parameters and five common calibrated GOTM parameters were included in the SA (Appendix A).The parameter space was defined as ±25% of the lake-specific calibrated parameter values, although the parameter space was bound by model constraints and not allowed to cross model limits (e.g. a parameter describing a fraction ranging from 0 to 1 was bound within this range) (Appendix A).
All parameter values were assumed to be uniformly distributed within the parameter space.
parsac (Bolding and Bruggeman, 2020) is a Python package for sensitivity analysis and calibration, designed specifically for models that take a long time to run.It allows the work for a sensitivity analysis or calibration procedure to be distributed over multiple compute cores of a single workstation or high performance computing cluster, which may potentially considerably reduce the time taken to complete the analysis.parsac further emphasises traceability and reuse of model analyses, which is made possible by logging parameter values and summary statistics for each model simulation to a MySQL database (https://www.mysql.com,last accessed November 1, 2020).The package is model-independent as it supports different mechanisms to configure parameters (yaml files, Fortran namelists, attributes of Python objects) and its access to model outputs (e.g., reading spatiotemporally varying fields from NetCDF files or calling Python functions).These features make it particularly suitable for analysis of GOTM-FABM, which is configured through yaml files and writes outputs in NetCDF format.Sensitivity analysis in parsac is built on top of SALib (Sensitivity Analysis Library in Python), which supports a range of commonly used sensitivity analysis methods.
This study applied two global sensitivity methods, the density-based Borgonovo's δ-method and the widely used variance-based Sobol' method, to calculate first-and total-order indices (Fig. 1).Initially, a SA with Borgonovo's δ was applied to identify and separate influential from noninfluential parameters (often referred to as the parameter fixing setting or parameter screening).For this purpose, a "dummy parameter" having no influence on the model output, was included in the sampled parameter set.In theory, non-influential parameters in the Borgonovo's δ and Sobol' total-order indices have a value of zero, but due to numerical approximations in the SA algorithms, the non-influential parameters may obtain small non-zero sensitivity indices.The sensitivity index of the dummy parameter thus provided an indication of the approximation error of estimating sensitivity indices.All parameters with a sensitivity index above the dummy parameter's confidence interval were classified as influential.Similar to this dummy approach has been applied in other sensitivity studies with complex environmental models (e.g.Zadeh et al., 2017).The parameter space was sampled with Latin hypercube sampling (LHS) with a sample size of 100, generating in total 36,200 model simulations for each lake model, and confidence intervals were determined by bootstrap resampling of 100 samples and a confidence interval level of 0.95 (Herman and Usher, 2017).
For further in-depth analysis of the shallowest and the deepest lake, the most influential parameters from Lake Hinge and Lake Ravn were selected based on the identified influential parameters from the screening analysis for the output variables: total chlorophyll a (chl-a), total phosphorus (TP) and total nitrogen (TN).To increase the strength of the analysis, a sample size of 250 was applied for both the Borgonovo's and Sobol' methods, generated with LHS and Saltelli's sampling scheme (Saltelli, 2002), respectively.To avoid numerical instabilities and consequently model crashes when running simulations with the generated parameter sets, the parameter sampling range for the in-depth analysis had to be constrained for Lake Ravn to ±15% of lake-specific calibrated parameter values.In the Lake Hinge in-depth SA, both ±25% and ±15% of lake-specific calibrated parameter values were applied to also compare in-depth and screening results.All results were analysed with confidence intervals determined by bootstrap resampling of 100 samples and a confidence interval level of 0.95 (Herman and Usher, 2017).

The density-based Borgonovo's δ-sensitivity measure
The density-based Borgonovo's method (Borgonovo, 2007;Plischke et al., 2013) is a global model and moment-independent sensitivity measure that considers the entire model output distribution in order to quantify the relative influence of parameters on the model output.The δ-sensitivity measure is based on estimating the change between the probability distribution function (PDF) of the model output and the PDF of the model output conditional upon fixing a given input parameter.
Let the model be represented by a function where Y is scalar model output, and X = (x 1 , …, x n ) is a set of n independent model parameters (in this study the model parameters).The separation between the unconditional probability distribution function of the model output (f Y (y)) and the probability density conditional upon fixing input parameter X i (f Y|Xi (y)) can be obtained by calculating The expected separation is then and the sensitivity measure δ i for parameter X i can be defined as The δ represents the normalised expected shift in the distribution of Y provoked by X i .The δ-sensitivity measure describes direct and indirect effects of the parameter and possesses several important properties: 1) δ lies between 0 and 1, 2) δ is zero when Y is independent of X i , and 3) δ equals unity when uncertainty in all model inputs is resolved (δ 1,2, …,n = 1) (Borgonovo, 2017).To estimate δ for model parameters, we applied the method laid out in Plischke et al. (2013) as included in the global sensitivity analysis software SALib (Herman and Usher, 2017).This method also includes calculating Sobol' first-order indices from LHS, so Sobol' first-order index is also applied to assess the screening results (e. g.Fig. 4).Model parameters were sampled with LHS, assuming a uniform parameter distribution with a sampling size of n = 100 and n = 250 for the screening analysis and in-depth analysis, respectively.One model output, Green algae DW, in the Lake Ravn model was log-transformed before Borgonovo's δ was estimated to handle unreasonably high δ-sensitivity measures.

The variance-based Sobol's indices
Sobol's sensitivity analysis method (Sobol, 2001;Saltelli, 2002;Saltelli et al., 2010) is a global and model independent sensitivity analysis that is based on variance decomposition.It can handle non-linear and non-monotonic functions and models.The variance-based Sobol's indices represent the fraction of the variance in the output that can be attributed to a given input alone (first-order: S i ) or with consideration of the input interaction with all other model parameters included in the analysis (total-order: S Ti ).
Let the model be represented by a function where Y is scalar model output, and X = (x 1 , …, x n ) is a set of n independent model parameters.The total unconditional variance of the model output V(Y) can be decomposed as where V i is the partial variance of X i on Y and is given by , while V ij is the impact of X i and X j on the total variance minus their first-order effects.Using these partial variances, the firstand total-order Sobol' sensitivity indices, S i and S Ti , can be defined in relation to the total variance (Borgonovo, 2017;Saltelli et al., 2008) where i = 1, …, n: Fig. 1.Illustrated workflow of the sensitivity analysis for each of the three lake model studies.The screening analysis included all FABM-PCLake parameters and five GOTM parameters in the density-based method, Borgonovo's δ-measure, sampled with Latin hybercube sampling (LHS).A dummy parameter was included to set a threshold in order to separate influential from non-influential parameters.If influential parameters were influential on two of the three model outputs: summer mean total phosphorous, total nitrogen and chlorophyll-a concentrations, these were included in an in-depth analysis applying both density-and variance-based methods, Borgonovo's δ-measure and Sobol' total-order indices.
To sample the parameter ranges and generate model inputs, this study applied the Saltelli extension of the Sobol sequence to estimate the first-and total-order indices at a cost of N = n(2p +2) evaluations, where N is the number of model simulations, n is the number of samples to generate, and p is the number of model parameters (Saltelli, 2002).Sobol' first-order sensitivity index is widely used and is a robust measure of functional dependence (Plischke et al., 2013;Saltelli et al., 2008).A null value of the first-order sensitivity measure does not imply that the model output is independent of the input parameter X i , as interactions are not captured (Plischke et al., 2013;Saltelli et al., 2008).

The lake ecosystem model: GOTM-FABM-PCLake
We applied a lake-adapted version of the one-dimensional water column model, GOTM, coupled with the lake ecosystem model FABM-PCLake.The hydrodynamic and ecological lake model, GOTM-FABM-PCLake, describes conceptually the most important hydrodynamic and ecological processes in temperate lake ecosystems.Based on lakespecific meteorological forcing, in-and outflows and a hypsograph (i.e. the relation between depth and the corresponding horizontal area), GOTM simulates the vertical transport of momentum, salt and heat through the water column.This study utilised the k-epsilon model for calculating the turbulence closure and the first-order, positive-preserving and conservative Extended Modified Patankar ordinary differential equation scheme for source and sink dynamics (Bruggeman et al., 2007).
FABM-PCLake is a dynamic, complex lake ecosystem model, which simulates nutrient and food web dynamics within the water column and the top sediment in each water layer (Hu et al., 2016).The model is a redesign of the widely applied aquatic ecosystem model PCLake (Janse, 2005).The modelled lake ecosystem consists of variables describing oxygen dynamics, organic and inorganic fractions of nitrogen, phosphorus and silica, three different groups of phytoplankton, zooplankton, zoobenthos, plankti-and benthivorous fish (juveniles and adults), piscivorous fish as well submerged vegetation (Fig. 2).A full list of model parameters and state-variables is available in Hu et al. (2016).

The lake model cases
A GOTM-FABM-PCLake model complex was set up, calibrated and validated for three Danish lakes (Lake Hinge, Lake Bryrup and Lake Ravn) during a lake model pilot project for the Danish Environmental Protection Agency.For calibrated parameter values, see Appendix A. These lakes range from shallow to deep in a Danish context.Lake Hinge has a mean and maximum depth of 1.2 m and 2.6 m, respectively, a surface area of 0.91 km 2 and a short hydraulic residence time of 16-32 days.The water column in Lake Hinge rarely stratifies for periods longer than a week during summer.The catchment area (approximately 55 km 2 ) comprises mainly agriculture (93%) (Nielsen et al., 2000), as reflected by the high external TP and TN loading of 2.1-3.7 Mg/yr and 62-121 Mg/yr, respectively, in 2001-2005.The lake is eutrophic with mean summer TP and chl-a concentrations of approx.0.1-0.2mg/L and 60-150 μg/L, respectively.A GOTM-FABM-PCLake model complex was set up, calibrated and validated on Lake Hinge for the period 1993-2007 by Andersen et al. (2020).
Lake Bryrup has a mean and maximum depth of 4.6 m and 9 m, a surface area of 0.37 km 2 and its hydraulic residence time varies between 2 and 3 months.The lake is polymictic and stratifies temporarily during summer.The catchment area is approximately 50 km 2 with 60% of the area being used for agriculture (Nielsen et al., 2000).The external loading of TP and TN to the lake was 0.7-1.2Mg/yr and 51-80 Mg/yr, respectively, in 2001-2005.The lake frequently experienced cyanobacteria summer blooms and had a mean summer TP and chl-a concentrations of approx.0.05-0.1 mg/L and 25-50 μg/L, respectively.A GOTM-FABM-PCLake model complex was set up, calibrated and validated on Lake Bryrup for the period 1996-2005 by Chen et al. (2019).
Lake Ravn has a mean and maximum depth of 15 m and 33 m, respectively, a surface area of 1.8 km 2 and a hydraulic residence time of approx.1.8 years.Lake Ravn is dimictic and the water column stratifies throughout summer until autumn (May to October-November).The lake catchment (57.2 km 2 ) is dominated by agriculture (70%) (Nielsen et al., 2000).In 2001-2005, the external nutrient loading of TP and TN to the lake was 0.8 Mg TP/yr and 47-108 Mg TN/yr.The lake is classified as mesotrophic with mean summer TP and chl-a concentrations of approx.15-23 μg/L and 9-12 μg/L, respectively.Based on a previous lake model by (Trolle et al., 2008a), a GOTM-FABM-PCLake model complex was set up, calibrated and validated on Lake Ravn for the period 1996-2005, with improvements in model performance (for calibrated parameter values, see Appendix A).
Meteorological data for atmospheric forcing of all three lake models were derived from the European ECMWF Interim dataset (Dee et al., 2011).Monthly averages of water inflow (m 3 /s) and nutrient concentrations (mg/L) were used as boundary conditions based on bimonthly measurements of inflow, inorganic and total nitrogen and phosphorous concentrations from the Danish monitoring program (NOVANA) for lakes (Johansson et al., 2015) combined with estimates from the DK-QNP model for the ungauged catchment area (Thodsen et al., 2016).The particulate organic fraction of the total phosphorus and nitrogen load was estimated for each individual lake as were the initialisation values for the sediment nutrient pool.

Screening by Borgonovo's δ-sensitivity measure
A total of 108,300 model simulations were executed during the screening step and analysed with the density-based Borgonovo's method for the three lake model studies.The screening step categorised 46, 23 and 21 parameters to be influential on two of the three output variables total chlorophyll a, TN and TP for each lake model for Lake Hinge, Lake Bryrup and Lake Ravn, respectively (Table 2, for full table see Appendix B).In total, 60 unique parameters were identified for further, in-depth analysis.

Parameter sensitivity across lake types
In the screening step, aggregated sensitivities of model parameters to model module level were overall similar for lake model outputs across the three different lake model studies with both the density-based δ-sensitivity measures (Fig. 3) and the variance-based Sobol' firstorder indices (Fig. 4).However, the models simulating the two deeper lakes, Lake Bryrup and Lake Ravn, were generally more sensitive to parameters related to mixing of the water column (e.g.wind_factor) for summer mean water temperature, nitrogen concentrations and phytoand zooplankton biomasses.In contrast, parameters related to the macrophytes module and particle dynamics (e.g.resuspension, kResu-sPhytMax or kTurbFish, and phytoplankton growth, cMuMaxGren) were more influential on determining model outputs in Lake Hinge compared with the two deeper lakes.However, as for TP both Lake Hinge and Lake Bryrup were more sensitive to parameters relating the benthic-pelagic coupling (e.g.resuspension) than the deepest lake, Lake Ravn.Thus, Lake Bryrup shows an intermediate response between the shallowest and deepest lake, depending on the model output.
The density-based sensitivity screening showed overall similar sensitivity patterns, though there were also clear differences in parameter sensitivity for some model outputs between the lakes.For all three lake models, a scale factor to the forcing wind speed (wind_factor) in the hydrodynamic model was the most influential parameter in determining the summer mean surface water temperature (Figs. 3 and 4).For the two deeper lakes (Lake Bryrup and Lake Ravn), where stratification occurs every summer, almost all of the variance of the water temperature could be ascribed to this particular scale factor (S1 = 0.95 ± 0.001 and 0.96 ± 0.0009, respectively).In the lake model for shallow Lake Hinge, the scale factor was also the most influential parameter (Hinge model δ = 0.224 ± 0003 and S1 = 0.125 ± 0.007); however, many parameters related to plankton and particle dynamics, such as growth rates of all three phytoplankton groups, zooplankton grazing and parameters related to particle resuspension by fish, together contributed with more variability to water temperature (Fig. 5).
For inorganic and total nitrogen concentrations, the Lake Bryrup and Lake Ravn models varied in influential parameters and modules compared with Lake Hinge.Phytoplankton growth and N-uptake parameters contributed most to the aggregated variability per module for outputs of summer mean nitrate (Lake Hinge and Lake Bryrup) and

Table 2
Ranking of the five most influential model parameters on total chlorophyll-a (chl-a), total phosphorous (TP) and total nitrogen (TN) from the sensitivity analysis screening step for each lake model.Model modules phytoplankton_water and abiotic_sediment are abbreviated to phyt_water and abiotic_sed.For a description of parameters, see Appendix A. ammonium concentrations (all three lakes).For the stratified lakes (Lake Bryrup and Lake Ravn), the nitrification rate constant (cNitrW) was the single most influential parameter explaining approx.30% of the variance in summer mean ammonium concentrations, when disregarding parameter interactions (Appendix B).Several zooplankton parameters (especially cFiltMax) influenced the summer mean nitrate concentrations in the Lake Hinge and Lake Ravn models (Fig. 3).Summer mean TN concentrations were mostly driven by parameters controlling particle concentrations (e.g.phytoplankton growth rates and zooplankton grazing) in Lake Hinge, parameters controlling phytoplankton growth and N mineralisation-related parameters as well as wind_factor in Lake Bryrup, and parameters related zooplankton grazing and to nitrogen release from the bottom sediments in Lake Ravn (Fig. 3).Interestingly, both the impacts of wind_factor and sediment parameters increased with the depth of the lake model.In the model for the deepest, stratified lake, wind_factor contributed with approx.28% of variance in the TN summer mean (without interaction effects).
For summer mean phosphate and total phosphorous concentrations, the parameter sensitivity response was overall similar for the three lake models.In all models, parameters in the zooplankton module were of primary importance and those in the phytoplankton module of secondary importance (except the Lake Bryrup TP output) in determining P concentrations.However, notable differences occurred.For Lake Hinge and Lake Ravn, parameters determining zooplankton grazing ranked highest in sensitivity of TP concentrations, and parameters describing P sediment mineralisation processes (e.g.kPMinPOMS and cTheta-MinPOMS) were also classified as influential.In the Lake Bryrup model, and also secondary for Lake Hinge, sedimentation and resuspension parameters were the most important in determining summer mean TP concentrations.
Phyto-and zooplankton outputs in general showed the same sensitivity response across all three lake models; however, some differences related to lake depth and trophic state occurred.Zooplankton grazing parameters (cFiltMax and hFilt) were the most dominant in determining summer mean phyto-and zooplankton biomasses for all lakes.For phytoplankton biomasses, maximum growth rates for each specific phytoplankton (cMuMaxPhyt) and other phytoplankton-related parameters (e.g.temperature modifiers) were also clearly influential (Fig. 5).

Fig. 3.
Borgonovo's δ-sensitivity measure for influential parameters aggregated to module level and scaled to one for each output variable for each lake model from the screening step with LHS.For a description of output variables see Table 1.(DW, dry-weight).
For most phytoplankton groups, the maximum growth rate of other phytoplankton groups was also influential; for instance, in the Lake Ravn model the most influential parameter on diatom DW was maximum growth rate of cyanobacteria (cMuMaxBlue).Silicate to DW ratios (cSiDDiat in both water and sediment phytoplankton modules) were the most influential parameters for both summer mean diatom and green algae biomasses in the Lake Bryrup model.In contrast to the deeper lakes, Lake Hinge had several influential parameters on phytoand zooplankton dynamics related to particle dynamics (phytoplankton growth rate, sedimentation and resuspension), highlighting the interaction between the sediment and the pelagic domains in driving plankton dynamics in this, the most shallow, lake.
Macrophyte growth rate (cMuMaxVeg) was the most influential parameter in determining the macrophyte biomass in the Lake Ravn and Lake Hinge models but was less important in the Lake Bryrup model.Here, the migration rate (cDVegIn) was the most influential parameter as a low calibrated value of the macrophyte growth rate prevented establishment of macrophytes.The ability of the models to simulate the light environment was also influential on macrophyte biomass as particle dynamic parameters (e.g.zooplankton grazing) and the non-visible fraction of shortwave radiation (A; GOTM parameter) were influential in turbid Lake Hinge and clear-water Lake Ravn, respectively.In Lake Hinge, the maximum growth rate of macrophytes was ranked in top 10 most influential parameters on summer mean ammonium, nitrate, TN, phosphorus, DO and macrophytes biomass, whereas maximum growth rate was influential in the Lake Bryrup and Lake Ravn models only for macrophytes and piscivorous biomass (Fig. 5).
The three fish groups were overall determined by parameters related to their food availability and predation from piscivorous fish on zooplankti-benthivorous fish.For zooplanktivorous fish (juvenile fish), the most influential parameters were related to phyto-and zooplankton, whereas for zoobenthivorous fish (adult fish), the single most influential parameter was zoobenthos assimilation efficiency (kDAssBent).The inclusion of a positive effect on piscivorous assimilation with increased macrophytes biomass in the FABM-PCLake model was illustrated in the Lake Ravn model where macrophyte growth rate (cMuMaxVeg, δ = 0.17 ± 0.005) was the most influential parameter on the piscivorous fish biomass.A similar response was not observed for the Lake Hinge.For a description of output variables see Table 1.Surprisingly, the most influential parameter in determining summer mean piscivorous fish biomass was the migration factor (cDPiscIn), and no predation parameters were found to be influential for Lake Bryrup.

In-depth analysis with Borgonovo's δ and Sobol' total-order indices
The in-depth analysis corroborated the initial findings of the screening density-based analysis for the two analysed models of shallow Lake Hinge and deep Lake Ravn.The relative ranking of the most influential parameters was consistent across the δ and Sobol' total-order indices, though the variance-based method gave consistently larger values when subtracting the δ measure of the dummy parameter from the estimated δ measures (for summer mean chlorophyll a concentrations see Fig. 6 and Appendix B).Zooplankton maximum filtration rate (cFiltMax) was the most influential parameter in determining summer mean TP and chlorophyll a concentrations for both Lake Hinge and Lake Ravn.Surprisingly, maximum growth rate of green algae was the most influential parameter for summer mean TN concentrations in Lake Hinge.The green algae group contributed with a minor fraction of the summer total phytoplankton biomass in Lake Hinge, and therefore calibration and validation of green algae-related model parameters were of less focus (e.g.N:DW ratios) (Andersen et al., 2020).The scale factor for wind speed, wind_factor, was the most influential parameter in determining summer mean TN concentrations.Other top ranked influential parameters on TP, TN and chlorophyll a concentrations were zooplankton preference factor for particulate organic matter (cPref-POM), zooplankton half-saturation constant (hFilt), maximum growth rate of diatoms (cMuMaxDiat) and temperature optimum for  1 and Appendix B, respectively.zooplankton (cTmOptZoo).

Parameter influence can differ between lakes
We undertook a global SA with both density-and variance-based sensitivity methods to analyse the sensitivity of model outputs in the 1D hydrodynamic lake ecosystem model GOTM-FABM-PCLake for three different lakes.Overall, the sensitivity responses of most summer mean model outputs were similar across the three lakes; however, some differences emerged, which can be attributed to the type of the lakes and specifically their differing morphology, ranging from shallow to deep.
The SA revealed increased influence of parameters describing the benthic-pelagic coupling for the shallowest modelled lake, Lake Hinge (Fig. 5).Resuspension and sedimentation parameters (wind_factor and all auxiliary parameters in Fig. 5) contributed with more variability to the summer mean temperature and nitrogen concentrations in Lake Hinge compared with the two deeper lakes.This corroborates the mechanistic description of physical processes in GOTM as less wind energy is required to break the water column stability and as wind events more easily resuspend particulates from the sediment in shallower than in deeper lakes (Kalff, 2002).The influence of resuspension as a main process in regulating particle dynamics with consequences for water turbidity, nutrient dynamics and biotic communities was also shown to be of key importance in another sensitivity study of a shallow lake model (Janse et al., 2010).Lake Hinge is eutrophic, so the amount of suspended particulates (i.e.high concentration of phytoplankton) also impacted light penetration and hence the heat captured in the water column, which likely also increased the influence of resuspension on the modelled temperature output.
The SA of the models of the deep Lake Ravn and to a lesser degree also of the less deep Lake Bryrup emphasised the importance of the wind_factor for temperature and its influence on the variability in nitrate and total nitrogen concentrations.The scale factor for wind modifies the air-water interactions that impact the mixing and water temperature in the epilimnion as well as the depth of the simulated thermocline.The parameter was calibrated for all three lakes to account for local variations that is not taken into account in the wind forcings, for instance a forest surrounding Lake Ravn (Trolle et al., 2008a).Both lake models correctly simulated summer build-up of phosphorus and ammonium concentrations in the hypolimnion (Chen et al., 2019;Trolle et al., 2008b).Variability in summer nitrate concentrations in surface waters could partly be explained by changes in thermocline depth, influencing hypolimnion ammonium concentrations and the transport of ammonium to the epilimnion where it is further oxidised to nitrate.The importance of the wind scale factor and the influence of the light parameters from the hydrodynamic GOTM model on many of the outputs of the FABM-PCLake model clearly illustrate the necessity of including the entire model complex and all parameters when performing a global SA to properly capture model sensitivity.
In general, the most influential parameters were closely related, but not limited, to their derivative output variable and many parameters contributed with variability in the modelled output in all three lake models.For example, summer mean phytoplankton, macrophyte and zooplankton biomasses were most sensitive to the maximum phytoplankton growth rate, the maximum macrophyte growth rate and the maximum zooplankton filtration rate, respectively.The same general conclusions have been drawn in other lake ecosystem model sensitivity studies (Janse et al., 2010;Makler-Pick et al., 2011), which stresses the Fig. 6.Comparison of screening (upper plots) and in-depth (lower plots) sensitivity indices on simulated summer mean chlorophyll-a concentrations for influential parameters for the Lake Hinge model with 25% parameter ranges.Borgonovo's δ-sensitivity measure with and without cut-off threshold from dummy parameter (upper left and upper right plots, respectively) from the screening step and Borgonovo's δ-sensitivity measure with cut-off threshold (lower left plot) and Sobol' totalorder indices (lower right plot) from the in-depth analysis.95% confidence intervals were estimated with 100 resample size.For corresponding parameter number and parameter name, see Appendix B.
importance of including all, not just a subset, of parameters in a SA.
In the in-depth analysis, zooplankton parameters such as maximum grazing rate (cFiltMax) and half-saturation of food (hFilt) were found to be most influential on chlorophyll a, TN and TP concentrations.Previous sensitivity studies of shallow lakes with the 0-dimensional (0D) PCLake model also ranked zooplankton grazing parameters as being very influential on simulated water quality variables (Janse et al., 2010;Nielsen et al., 2014).Both 1D FABM-PCLake and 0D PCLake include only one zooplankton group in the food web.In comparison, a hydrodynamic-ecological lake model, DYRESM-CAEDYM, applied to Lake Kinneret, Israel, included three groups of zooplankton, and a global sensitivity study of this model did not rank zooplankton grazing parameters as most influential, although the messy feeding and grazing rate of zooplankton was overall influential (Makler-Pick et al., 2011).The secondary importance of zooplankton parameters compared with nutrient-related processes may be explained by the fact that Lake Kinneret is bottom-up controlled (Makler-Pick et al., 2011) as is typical for warm lakes due to the overall higher fish predation on zooplankton suppressing the abundance and size (Jeppesen et al., 2020).Parameterisation of only one zooplankton group to describe an entire zooplankton community and its impact on the lake ecosystem is probably a too crude and insufficient assumption as most zooplankton communities vary with the season in both shallow and deep lakes across a trophic gradient (Jeppesen et al., 2011).Increasing the number of zooplankton groups would likely partition some of the sensitivity and increase the accuracy of the parameterisation.
The SA highlights the sensitivity of possible shifts in the simulated dominance of primary production between phytoplankton and macrophytes in shallow lakes by FABM-PCLake.Hence, the trophic state of lakes, which also differed between the three lakes in our study, will likely also shape the sensitivity of the models.Similar results were found in a local SA applying 0D PCLake to Lake Arreskov, Denmark (Nielsen et al., 2014).The maximum growth rate of the macrophytes (cMu-MaxVeg) was in the top 10 of the most influential parameters on summer mean chlorophyll a, TN, TP concentrations and macrophyte biomass (Nielsen et al., 2014).However, a SA conducted with the method of Morris (Janse et al., 2010) did not find cMuMaxVeg to be influential, not even on macrophyte biomass, although other macrophyte parameters, such as the overwintering fraction of submerged vegetation (fWinVeg) and macrophyte respiration rate (kDRespVeg), were shown to impact the variability of simulated TN and TP outputs (Janse et al., 2010).
This SA was confined to analyse parameter sensitivity within the conceptual structure of GOTM-FABM-PCLake and did not evaluate impacts of different model structures.For example, ice coverage was not included in the GOTM-FABM-PCLake version applied in this study.However, it is now possible to simulate ice coverage in GOTM (htt p://www.github.com/gotm-model,last accessed November 1, 2020), and it cannot be ruled out that this might have impacted both the simulated water quality and ecological state variables (Balayla et al., 2010) and, thereby, parameter sensitivity.With continued model development, this highlights the need for more frequent application of SA in the modelling process.

Parameter sensitivity as a model evaluation tool
A common modelling procedure is to rely on previous sensitivity studies when deciding which parameters to target during calibration; and only rarely is a new SA conducted for a new model application to an individual lake.However, facilitated by parsac, the SA methods presented in this study can now, without much effort, be applied before any new lake model application undergoes calibration.Parameter ranges in an initial SA should be based on literature values, expert knowledge or similar model-set ups as there are no calibration results to depend on.
During model calibration, a SA could serve as an additional model evaluation tool and provide a system-level assessment (as presented in Hipsey et al. ( 2020)), under the assumption that the influential parameters or modules (when importance measures have been aggregated to module-level), identified by the SA of the models, reflect the most important processes in each of the lakes.Then, if a discrepancy was identified between empirical data and the outcomes of a SA in terms of the most influential parameters and processes, this would warrant an investigation into the specific process and/or entire model parameterisation and the cause of the discrepancy, and thereby help improve the understanding of the interactions in lake ecosystem.Potential discrepancies, could among other reasons be ascribed to wrongful parameter calibration, the chance of the calibration process converging on one among several plausible combinations of model parameter values (i.e.equifinality), parameter ranges applied in the SA, missing knowledge of processes in the ecosystem and/or conceptual errors in the lake models.If wrongful parameter values or equifinalty is the cause of the discrepancy, this could warrant a new round of calibration or changes to the model structure to better reflect the knowledge of influential processes in the lake ecosystem.For example, a decoupling of the simulated piscivorous fish from the rest of the food web was identified in Lake Bryrup, and a surprising impact of simulated green algae on TN concentrations was identified in Lake Hinge.The two examples are not in correspondence with current knowledge of the systems modelled and therefore warrant investigations into these parameterisations and maybe further calibrations and/or changes in the model structure.

Technical assessment of analysis approaches
With the newly developed software application parsac, we ran a global SA with both density-and variance-based sensitivity methods in a two-step approach (screening and in-depth) to fully analyse the hydrodynamic lake ecosystem model GOTM-FABM-PCLake.By utilising parallel processing in the execution of model simulations, parsac markedly decreased computing time, thereby removing one of the obstacles to applying continuous sensitivity studies in the model development and application process (Jakeman et al., 2006;Robson et al., 2008).For example, on a Microsoft Windows Server 2016 with four processors with 64 cores, executing in total 36,100 model simulations for Lake Bryrup took less than 5 h.parsac also benefits from the open source community as it is written in Python and utilises SALib (Sensitivity Analysis Library in Python) (Herman and Usher, 2017).

Variance-based sensitivity corroborates density-based parameter ranking
The two-step approach with a screening analysis followed by an indepth analysis, applied in this study, resembles other global SA approaches combining a screening and a global SA method as recommended by, for instance, Saltelli et al. (2008).By way of example, a sensitivity study of the biogeochemical lake model PCLake (Janse et al., 2010) utilised the Morris' method (Morris, 1991;Saltelli et al., 2008) and then the variance-based Fourier amplitude sensitivity test method (Saltelli et al., 2008), whereas the sensitivity study of the lake ecosystem model CAEDYM applied first recursive partitioning and regression trees (Breiman et al., 1984) and general linear model (Laird and Ware, 1982) analysis methods and, subsequently, a generalized boosted modelling analysis (Friedman, 1999).Here, we conducted an initial screening applying the density-based method Borgonovo's δ-sensitivity measure, with a lower sampling size, in combination with an in-depth analysis as this method met several criteria established for integrated models to be used in environmental management and decision-making (Ravalico et al., 2005): the method was capable of handling +300 parameters, taking into account higher order parameter interactions and model non-linearity, not requiring knowledge of parameter probability distributions and output interpretable results.Another important criterion was that the method should be able to handle eventual removal of parameter sets (i.e. for simulations that resulted in a model failure) as running a global SA without detecting model errors would be unusual (Saltelli et al., 2019).Borgonovo's δ-sensitivity measure relies on LHS, which allows for removal of parameter sets without breaking the sampling strategy (the removal of parameter sets was also tested, see Appendix B).Therefore, what is traditionally viewed as a screening step was in this study in reality a full global SA, and the following SA step was a test of the screening results with both a density-and variance-based method.
The density-and variance-based methods agreed in the ranking of the most influential parameters on simulated summer mean TN, TP and chlorophyll-a concentrations for the shallow (Lake Hinge) and deep (Lake Ravn) lakes subjected to in-depth analysis (Figure S6, Figure S7 and Figure S8).The agreement between the parameter ranking of the density-and variance-based methods increased the confidence in which parameters are important (Borgonovo et al., 2017).Thus, applying several SA methods following the initial density-based analysis will function as a test and possible corroboration of the density-based analysis.Instead of applying other sampling strategies and corresponding sensitivity methods, several sensitivity methods can be based on LHS as exemplified in the given-data approach by Borgonovo et al. (2017).For example, first-order Sobol', Borgonovo's δ, Kuiper metric and the DELSA methods were used to interrogate a hydrological model (Borgonovo et al., 2017).
Overall, the screening and in-depth analysis replicated the ranking of the most influential parameters between the two lakes.As the screening step in this study was a full global SA (though with a lower sampling size), and not a sensitivity screening method, the discrepancy between the screening step and the in-depth analysis step should be relatively low.An advantage of the screening analysis was the inclusion of all model parameters, which impacted the identification of most sensitive parameters in this study.

Impacts of parameter ranges on sensitivity
The definition of parameter ranges included in the SA potentially has a significant impact on the results and should therefore be considered carefully in all sensitivity studies (Wagener and Pianosi, 2019).Setting plausible ranges of allowed parameter values when sampling for the analysis may be a challenging task, however.Sensitivity studies of aquatic ecosystem models have based parameter ranges on both literature and expert judgment (e.g.Makler-Pick et al., 2011;Missaghi et al., 2013), percentage ranges from calibrated and/or default model values (e.g. Bruce et al., 2006;Elliott et al., 1999) or a combination of these (e.g.Omlin et al., 2001).To properly capture the targeted lake's specific ecosystem dynamics and simulate different lake types, this study applied ±25% and ±15% ranges on calibrated and validated parameter values.For the Lake Hinge model, total-order Sobol' and Borgonovo's δ in-depth analysis agreed in the ranking of the most influential parameters on summer TN, TP and chlorophyll a concentrations with 25% and ±15% ranges on calibrated parameter values.In contrast, the Lake Ravn model was impacted by parameter range size as only ±15% ranges, and not ±25% ranges, sampled with Saltelli sampling allowed error-free model simulations, but this does not necessarily mean that the parameter sensitivity would be different if model execution with 25% ranges were possible.We also searched for patterns in the parameter sets of model crashes with a decision tree analysis (DecisionTreeClassifier in scikit-learn version 0.23.1 (Pedregosa et al., 2011), results not shown); however no clear patterns emerged, so removal of model crashes and their parameter sets likely did not impact the results.
As GOTM-FABM-PCLake is a complex model characterized with the equifinality issue (Beven, 2006), the SA results could potentially also be impacted by this issue as the calibrated parameter values served as basis for parameter ranges in the SA, even though the parameter range size of ±25% and ±15% was applied.The equifinality issue re-emphasises the need to not rely on a single or few SA studies, but to apply SA actively in the modelling process.

Conclusions
Parallel Sensitivity and Auto-Calibration (parsac) with density-and variance-based sensitivity methods offers a powerful tool to aquatic ecosystem modellers investigating numerical models by applying global SA.Applying the density-based global SA, Borgonovo's method, as a screening step proved efficient in ranking the most influential parameters, and the ranking was corroborated by the variance-based SA, Sobol' method.When performing a global SA with lake ecosystem models, all model parameters and the entire model complex should initially be included in the analysis, and SA should preferably be performed for each individual study case as lake morphology and lake type can shape the parameter sensitivity of lake ecosystem models.For total phosphorus, we found that parameters related to the benthic-pelagic coupling were particularly important for the shallower lakes and less so for the deepest lake.For total nitrogen, and to a lesser degree for nitrate, we found that phytoplankton-related parameters were the most sensitive for the shallower lakes, whereas parameters related to the benthic-pelagic coupling were more important for the deepest lake.For chlorophyll a and the biomass of all phytoplankton groups, we generally found that phytoplankton and zooplankton parameters were most important, the latter to a lesser extent in the shallowest lake where resuspension also impacted output variability.We infer that lake morphology shapes the parameter sensitivity of lake ecosystem models and should be factored in when experimenting with these models.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 4 .
Fig. 4. Sobol' first-order indices for influential parameters aggregated to module level for each output variable for each lake model from the screening step with LHS.For a description of output variables see Table1.

Fig. 5 .
Fig. 5. Pattern plot of each lake model's influential parameter's Borgonovo's δ-sensitivity measure for all included output variables from the list of 60 most important overall parameters identified in the screening step.For a description of output variables and model parameters seeTable 1 and Appendix B, respectively.
T.K.Andersen et al.