1 INTRODUCTION

Marine ecosystemdynamics has become increasingly importantin the field of oceanographic research. Despite its complexity, our understanding of ecosystem dynamicprocesses is expandingfollowing improvementsin observational and experimental capabilities (Friedrichs, 2002; Tang et al., 2005; Shi et al., 2012; Zheng et al., 2012) and developments from single disciplines, such as physical oceanography, marine chemistry, and marine biology, and interdisciplinary fields such as marineecology (Rose et al., 2010) . Marine ecosystem dynamic models (MEDMs) are important tools with which the dynamic processes of marine ecosystems may be described quantitatively. Recently, MEDMs have been applied extensively to studies of the structureand function of marine ecosystems, to forecast theirhealth and to assess their ability to respond to global climate changes (Chen et al., 1999; Friedrichs et al., 2007) . MEDMs are groups of mathematic functionswith multiple parameters and state variables, which have increasing nonlinearity and complexity (Friedrichs, 2002; Tang et al., 2005; Fulton, 2010; Shi et al., 2011, 2012; Zheng et al., 2012) . The common structureof a MEDM encompasses the lower trophic levels, i.e., usually, the classical model of nutrient-phytoplankton-zooplankton. Thus, the components of organic detritus, bacteria, dissolved organic material, and aquaculture organisms such as algae, shellfish, and fish are added according to modelers’ specific interests (Rose et al., 2010) . Models of lower trophic levels are generally dividedinto components, which comprisemultiple functional groups of different sizes and species; however, there is an increasing requirement for model applications to consider higher trophic levels, particularly anthropogeniccomponents (Fulton, 2010; Rose et al., 2010; Morris et al., 2014; Xiao and Friedrichs, 2014a) . Thus, it is necessary for MEDMs to exhibit good model skill in forecasting dynamic changes of marine ecosystems and to provide a scientific basis for marine resources management, especially ecosystem-based marine management (Stow et al., 2009) . Therefore, the improvement and assessment of the skill of MEDMs is an important issue.

A general definition of the term “model skill”is presented based on the discussion of skill assessment for coupled biological/physical models of marine systems, which was sponsoredby the Journal of Marine Systemsin 2009. Therein, model skill is specified as the fidelityof model behaviorto natural truth and assessed basedon human judgment. Because natural truth is generally approximated by observational results, model skill can be regarded as how well the model fits the observational data. Thus, the suitability of model skill becomes the basis on which tojudge the reliability of the application of the model. It is imperative that model skill together with the reliability and robustness of its application be optimized and improved (Stow et al., 2009) . Moreover, quantitative assessment of model skill could be used to direct its optimization andapplication.

The primary factorsaffecting model skill are the model structure, parameters, solution, and inputs (Gregg et al., 2009) . Among these, model structureis distinct from the others becauseit principally refersto the biological process, coupling mechanism with other processessuch as physical and chemical processes, spatial dimensionality, and model complexity. Because the dynamical processof the marine ecosystem in the model is based on hydrodynamics, its coupling mechanismand spatial dimensionality depend upon how the dynamical process is affectedby hydrodynamic conditions such as horizontal advectionand vertical mixing.Model complexitygenerally refers to the numberof state variables and it depends on the type of modeled ecosystems (Arhonditsis and Brett, 2004) . Following the development of physical processresearch, the currently acceptedpractice is that the biologicalmode of low complexity is coupled with the integrated physical mode. Generally, it is believedthat by increasing the complexity of the modelstructure, such as couplingboth physical and chemical processes, using three-dimensional modes, and increasing the number of state variables, the potential of the model to reproduce a complex naturalsystem accuratelyis increased. However, it has been proven that modelskill is not reliant on the type of modeled ecosystem and that an increase in the complexity of model structure does not necessarily lead to a systematic improvement in model skill (Matear, 1995; Hurtt and Armstrong, 1996, 1999; Arhonditsis and Brett, 2004; Flynn, 2005; Friedrichs et al., 2007, 2009; Xiao and Friedrichs, 2014a) . The lack of full knowledge of a complexmarine ecosystem would lead to insuffcient support of a model with relatively high complexity, i.e., the additional state variables and parameters would inevitably introducesignificant uncertainties. Accordingly, there has to be trade-offs, principally between modelstructure and model skill (Ward et al., 2013; Xiao and Friedrichs, 2014a) .Thus, skill assessment can be used to assist modelers in theirselection of modelstructure. The othertwo factors of model solution and inputs mainly depend on specific numerical algorithms and observational data, respectively. Nevertheless, the exact values of model parameters are diffcult to determine or measure because they vary with differentmodels and study areas, which means there is no generallaw to follow (Friedrichs, 2002) . Ultimately, it can be concluded that the impact factorsof model structure, solution, and inputs dependon the development level of theoretical research and observational technology, while the selection of model parametersremains relatively flexible.For this reason, an improvement in the skill of MEDMs is generally dependenton parameter optimization (PO) . There is wide recognition that modelparameters can be optimized mathematically, mainly because a model is always a conceptual expression of natural processrather than an actualcopy. That means, for example, multiple classes and species of different sizes may be classified into a functional group, such as phytoplankton and zooplankton, based on which the evolution of different state variables and their intercommunications can be described by empirical functions (Hood et al., 2006; Peng et al., 2012) . However, the parameters of these empirical functionsare always obtainedthrough laboratory measurements on individual species under controlled conditions and therefore, it is diffcult to use them to describethe overall biological characteristics of a marine ecosystem (Pahlow and Oschlies, 2009) . Moreover, experimental data are limited in time and space, such that the value of the same parametercan vary considerably between different literatures; thus, researchers cannot obtain definite parametervalues but only a rational range. Previously, model parameters have been estimated through manual adjustment according to both empirical accumulation and subjective judgment, which is tedious work and makes it diffcult to achieve satisfying results.Therefore, to improvemodel skill, automatic optimization methods for improvingthe goodness-of-fit of simulations to observations are often used to adjust model parameters (Gregg et al., 2009; Ward et al., 2010) .

The PO problem related to MEDMs is often nonlinear, multi-objective, and multi-extremum, which presentsmany diffculties in optimization. First, becausethe large numbers of parameters related to nonlinearmodels need tremendous solution space for the optimization problem, the amount of observational data and the search abilityof the optimization algorithm must obey high requirements. Second, the optimization problem refers to multiple objective functionsbecause of the multiple state variables. Third, with the increasein model complexity, the analyticproperties of the objective functions requiredby common optimization methods are unavailable. For example, even with high search abilities, gradient methods remain unsuitable for complex modelsbecause of the demand for the gradient information of objective functions. Finally, the nonlinearity and complexity of MEDMs make the relevant goal function a multimodalfunction, which causes the optimization problem to become a multi- extremum problem.Nevertheless, for the majority of optimization methods, such as the simplex method, because the search processhas no strong randomness and its search ability depends greatly on the starting point, it is easily trapped in local extrema, which can lead to premature convergence and thus, the missing of a global optimum (Gregg et al., 2009) . These cases would become increasingly severe with the increase of the solution space of the optimization problem.

To facilitate skill improvement of MEDMs, efforts have been made to explore basic procedures, as outlined in the frameworkshown in Fig. 1. First, to solve the diffcultiesmentioned above, two key problems must be addressed regarding the parameterization of MEDMs. The first is to identify those parameters that have the greatest effect on the model results, and the other is to establishhow to optimize them to improve the goodness-of-fit of the simulationsto the observations. The former can be determined through parameter sensitivity analysis (SA) and the latter relies on the effectiveselection of an optimization algorithm. Therefore, the PO based on SA (generally called model calibration) is not only the primarystep but is also a challenge in modelling (Fennel et al., 2001; Solidoro et al., 2003) . Subsequent model validation is also an essential procedure (Gregg et al., 2009) . This step must be performedto examine the skills of the calibratedmodels to provide a scientific judgment on the effciencyof the optimized models. This is mainlybecause the optimized models are likely to absorb those errors of the model parameters introduced from physical models or other unknown processesduring the procedureof PO (Gregg et al., 2009) .Additionally, model validation can assist in the evaluation of the skills of the PO algorithms.

Currently, even though researchon MEDMs is ongoing, the level of research on the improvement and assessment of skill variesand there is a lack of a comprehensive and in-depth review. This paper summarizes the methods and strategies used in the improvement and assessment of the skill of MEDMs; it may not cover every aspect, but the objectivewas to inspire furtherresearch and to attract greaterattention to model skill.The hope was to form a scientific and normative technical framework for the improvement of the skill of MEDMs, ratherthan to provide a manual.

2 PARAMETER OPTIMIZATION FOR MODEL CALIBRATION 2.1 Sensitivity analysisfor selection of control parameters

The huge solutionspace of the parameter optimization problemplaces high demandsin terms of the amount of observational data required and the search abilityof the optimization algorithms. However, because obtainingmarine observations requires considerable time and effort, in general, neither the observational data nor the optimization algorithms meet these requirements. Thus, decreasing the dimensions of the solutionspace is considered an effcient means to reduce the diffculty of the optimization problem.Broadly, even with the large number of parameters involvedin the MEDMs, only a limitednumber would have significant impact on the modelresults. Therefore, by analyzing the impacts of the parameters on the modelresult s, those parameters with the largest impactscould be selected as the control parameters for further optimization. Parameter SA is an effectivemethod with which to analyze the impacts of parameters on model results and it has been considered as a primarystep of PO in MEDMs (Fennel et al., 2001; Solidoro et al., 2003; Shi et al., 2012, 2014) . Furthermore, highlycorrelated parameterscannot be optimized effectively and simultaneously. Friedrichs et al. (2007) expounded that a change made to one parameterwould be counteracted by an associated change in another highly correlated one, resulting in multiple combinations of parameter values producing indistinguishable results during the procedure of PO. Thus, parameters with weak inter-correlations as well as large impactson the model results have been proven to comprise the best set of control parameters (Chu et al., 2007; Cossarini and Solidoro, 2008; Gibson and Spitz, 2011; Zheng et al., 2012; Morris et al., 2014) . Accordingto whether the interactions among parameters are taken into consideration, SA can be divided into two types:local sensitivity analysis (LSA) and global sensitivity analysis (GSA) .

2.1.1 Local sensitivity analysis (LSA)

LSA studies the impact of a single parameter with a small variation around a certainvalue. Herein, the sensitivity index is generallydefined as the change rate of the model result divided by the related change rate of the parameter. As a simple and operabletool, it has wide applications in the field of MEDMs (Solidoro et al., 2003; Shi et al., 2012, 2014) .However, it does have some limitations (Shi et al., 2012, 2014) . First, the LSA result relies on the initial value of the parameter and its variation. Second, artificiallytaking the variation as the fixed value might cause the values of some parameters to exceed their actual domains. Third, the impact of the interactions among the parameters on the analytic results cannot be teste d. Consequently, LSA might not be a useful approach for the selectionof control parameters. Nevertheless, some researchers have adopted the statistical approach to analyze LSA resultswhen trying to select an appropriate set of controlparameters (Solidoro et al., 2003; Shi et al., 2014) . Althoughits effectiveness has yet to be validated, this method could be used to examine the correlations amongthe parameters, and whencomputational capacity is limited, it could be regarded as a feasiblemethod.

2.1.2 Global sensitivity analysis (GSA)

GSA simultaneously tests the impacts of multiple parameters with variations in their actual ranges. As a global analysis, it has advantagesover LSA. On the one hand, parameters vary across their entire value domains, independent on their initial values. On the other hand, variationof the model results is caused by the joint variation of all the parameters and thus, the parameter sensitivity is global. Using GSA, those parameters with significant impact on the model results but no strong inter-correlations can be selected for the set of control parameters.

There are several methods for implementing GSA, as illustrated by Cariboni et al. (2007) . Among them, the Morris and Sobol methods are the two approaches used most frequently in the field of MEDMs (Chu-Agor et al., 2011; Morris et al., 2014) . The Morris method is a state-of-the-art screening approach. By measuring the relative impacts of the parameters on the model results, it provides a so-called qualitative ranking of the parameters (Cossarini and Solidoro, 2008; Zheng et al., 2012;Morris et al., 2014) . Conversely, a quantitative analysis of parameter sensitivity can be performed using the Sobol method (Sobol, 1993) as well as the extendedFourier amplitude sensitivity test (Saltelli et al., 1999) . It is a variance-base d approachand thus, it is capableof quantifying the amount of varianceexplained by arbitrary combinations of model parameters (Chu et al., 2007; Gibson and Spitz, 2011; Morris et al., 2014) .

Despite the limitations of GSA comparedwith LSA in terms of highercomputational cost, the advantage of its samplingtechniques makes it a more powerful and flexible tool with which to explorethe entire range of the variationof the parameters simultaneously (Sobol, 1993; Gibson and Spitz, 2011) . For example, the sampling technique designed by Morris (1991) and extended by Campolongo et al. (2007) is skillfuland computationally cheap. Additionally, Latin hypercubesampling is a stratified sampling scheme of the Monte Carlo simulation, which can avoid the samplingbeing aggregated within the range with high probability and thus, it is an effcient technique for the Sobol method (Iman, 2008) . To exploitthe Morris and Sobol methodsfully, common practiceis to first excludethose parameters with low sensitivity using the Morrismethod, and then to quantify the sensitivity of the remaining parameters using the Sobol method (Chu-Agor et al., 2011; Morris et al., 2014) .

2.2 PO methods for improvement of the goodness- of-fit

Model parameteroptimization is the use of existing observational data to inversemodel parameters; it alsocalled data assimilation (DA) in the fields of meteorology and oceanography. DA has been applied to optimize the initial conditions, open boundary conditions, and model parameters in marine numerical simulationsof tidal currentfields and sea temperatures. Since Ishizaka (1990) first assimilated Coastal Zone Color Scannerdata into a three-dimensional physical- biological model using data insertion, the use of DA has graduallyincreased in relationto MEDMs (Gregg et al., 2009) . According to the different sequencesof using observational data, DA is divided into inverse and sequential DA (Bouttier and Courtier, 1999; Gregg et al., 2009) . The inverse method minimizes the cost function to obtain an optimized parameter set. It includes the steepest descentmethod, conjugate gradient method, variational adjoint, and also meta- heuristic algorithms represented by the simulated annealing and genetic algorithm. The sequential method initializes the model by correcting its output at the time when observational data are available, and it contains the nudging method, optimal interpolation, and Kalman filter as well as its various transformations. Gregg et al. (2009) summarized the applications of DA to the field of marine biology. They highlighted that the inverse method has historically been the most popular for parameterization. The sequential method has been used less but its use has been growing in popularity since about 2000, and it is commonlyused in the estimation and prediction of state variables. Therefore, PO can be considered an effective means of DA.

The performance of optimization algorithmsis described by the convergence property and rate, which refer to the precision of the optimized models and the effciency of the algorithms, respectively. However, it is found that common optimization algorithms easily converge to their local optima, resulting in optimized models that are not suffciently precise. For example, from the designpoint of view, the commonvariational adjoint and Kalman filter algorithms are both local methods that are unableto resolve effectivelythe problem of multi-extrema optimization constrained by MEDMs, which means researchers have increasingly focusedon global optimization, i.e., algorithms with global search ability, such as the simulated annealingand genetic algorithm.

2.2.1 Common local PO methods

Currently, PO methods used for MEDMsare mainly local. This is because research on the PO of MEDMs has occurred latelyand local methodsare relativelymature, and the precision of some less nonlinear models can be improved greatly by local methods. Furthermore, global methods have not yet matured and their computational costs are high because of their low convergence rates. Accordingto the differentrequirements on the analytic properties of optimization problems, local methodsinclude two types: gradient-based numerical methods and direct searching methods (Tashkova et al., 2012) .Gradient- based numerical methods, which generally refer to inverse methodsexcluding meta-heuristic algorithms, are highly effcient but poor at tacklingproblems with discontinuous, unsmooth, polymorphic, or sick goal functions. Direct searching methods such as sequential methods are not very effcient, despite their low requirements on the analyticproperties of optimization problems. However, as a type of early and relatively mature algorithm, local methods have been widely used.

Because the searching effciency of global methods commonly depends on the local searchingcapability, the research of local methods contributes greatly to the development of global methods. Presently, the most widelyused algorithms are the variational adjoint and Kalman filter (Gregg et al., 2009) . Variational adjoint combinesobservational data and MEDMs using the variational method and adjoint function, adjustingthe control parameters according to the gradientsof the goal function. It was firstly used by Sasaki (1970) and was applied to the inverse problem of a tidal current model in open water s by Bennett and McIntosh (1982) . Since Lawson et al. (1995) adoptedthis method to deal with the parameters and initialfields of a simple marine predation-prey model, variational adjoint has been widely used in MEDMs (Zhao et al., 2005; Fan and Lv, 2009; Gregg et al., 2009; Ward et al., 2010) . However, these methods have limitations. For example, variational adjoint relies considerably on the initial values of the parameters, the spatiotemporal scale of the observational data, and it requiresgradients. The Kalman filter constantly updates the predictive field using new observational data as the subsequent initial field to improve the predictiveability of models; this can be regardedas an extension of optimal interpolation (Jones, 1965; Evensen, 1994; Gregg et al., 2009) . The difference of this method compared with others is that it can predict and analyze state variablesand errors simultaneously. The idea of the Kalman filter was proposed by Kalman (1960) and it was called the standard Kalman filter. Since Jones (1965) introduced it to meteorology, it has been applied widelyin the fields of both meteorology and oceanography. Following the improvement of model nonlinearity, the extended, simplified, and ensemble Kalman filters have graduallyemerged to improve the convergence and stabilityof the Kalman filter in relationto nonlinear systems.The ensemble Kalman filter was proposed by Epstein (1969) , based on the theory of random dynamicprediction, and it was applied to DA by Evensen (1994) . The ensemble Kalman filter uses ensemble statistics to estimate the error covariance and thus, is suitable for an optimization problem of a highly nonlinear model (Gregg et al., 2009) .

2.2.2 Exploring global PO methods

To deal with the optimization problem of multi- extrema MEDMs, global optimization approaches are used to improve the precisionof the model parameterization. According to whetherstate variables are considered as random variables, global optimizationapproaches can be classified as deterministic and stochasticmethods (Tashkova et al., 2012) . They have different requirements on the analytic properties of optimization problems. Based on the analytic propertyof the optimization problem, deterministic methodsproduce a deterministic finite or infinite point range with a sub-sequence converging to a global optimum, which theoretically ensuresthat the optimal solutionis global. However, the algorithm of these methods variesdepending on the specific analytic property of the optimization problem. Moreover, becausethe point range might be infinite, itis not guaranteed that any type of problem could be solved within a limited time. In addition, the computational cost increasesrapidly with the increase of the problem scale. Therefore, this approach is unsuitable for complex models despite its application in simple nonlinear dynamic system models in the fields of biology and chemistry (Athias et al., 2000; Esposito and Floudas, 2000; Miró et al., 2012) .

Theoretically, stochastic methods cannot guarantee to obtain a globaloptimal solution, and the degreeof deviation between the obtained optimalsolution and the global optimal solutionis unpredictable. However, based on a probabilistic search, this approachcan effciently locate a position adjacentto the global optimal solution with acceptable computational costs and obtain an approximate or quasi-global optimal solution (Athias et al., 2000; Tashkova et al., 2012) . Stochastic methods, which belong to direct searching methods, have no requirement (or no high requirement) to the analytic properties of the optimization problem and even no concern about whether its expression is explicit; thus, they play an irreplaceable role in the PO of complex models (Athias et al., 2000; Tashkova et al., 2012) . Stochastic methods comprise a series of meta-heuristic algorithms. They combine local searching with a stochastic approach, and they are generally based on the imitation of natural phenomena and social behaviors, includingevolutionary algorithms, simulated annealing, swarm intelligence, and tabu search (Table 1) . Currently, stochastic methods are widely used in models of watershed hydrology and water quality, but rarely in MEDMs, and they tend to be limitedto genetic algorithms (Schartau and Oschlies, 2003; Oschlies and Schartau, 2005; Weber et al., 2005; Ward et al., 2010) and simulated annealing (Matear, 1995; Hurtt and Armstrong, 1996, 1999; Kavanagh et al., 2004) .

Following the improvement of the nonlinearity of MEDMs and the increase of structural complexity, common meta-heuristic algorithms such as genetic algorithms and simulated annealing have gradually exposed the shortcomings of the poor capability of local searching, locally premature convergence, and non-satisfied rates of convergence.Therefore, related modification is necessary. For example, the standard genetic algorithm operates a large population, which means numerous searching points to avoidpremature convergence when dealing with a high-dimensional search space.Thus, it has high computational cost and time overheads, whereas the micro-genetic algorithm, which is one of the genetic algorithms used most widely in MEDMs, overcomessuch disadvantages using a small population and re-initialization scheme (Schartau and Oschlies, 2003; Oschlies and Schartau, 2005; Weber et al., 2005; Ward et al., 2010) .

Because there is no perfect theory with which to analyze the convergence and its rate of different algorithms, comparison of algorithm skill is commonly based on numericaltests. Because the optimization effciency of meta-heuristic algorithms relies considerably on the propertiesof the optimal solution itself and computational cost, different case studies result in different answers.Recent research in related fields has shown that the micro-genetic algorithm is superior to simulated annealing, and that evolution strategiesare superior to generic algorithms. For example, Athias et al. (2000) used a complexand nonlinear marin e biogeochemical model to simulate the cycle of the marine particle system and found that compared with simulated annealing, the micro- genetic algorithmshowed higher computational effciency, a faster convergence rate, and a more stable optimal solution. Moles et al. (2003) compared several global optimization methods in the parameter estimation of nonlinear dynamicbiochemical pathways, and concluded that evolution strategies might be the most competitive stochastic optimization method. Huret et al. (2007) also successfully applied evolution strategies to a phytoplankton production model of the Bay of Biscay. Except for the properties of the optimal solution itself, the optimization effciency relies considerably on computational capacity. According to the PO of a complex hydrologic model using algorithms such as the genetic algorithm, particle swarm optimization, differential evolution, and artificial immunesystem optimization, Zhang et al. (2009) suggested that the genetic algorithm should be chosen when computational capacity is available, but that particle swarm optimization should be preferred when computational capacityis limited.

Current research on stochastic globaloptimization methods is focusedon their realization, modification, and application although the relatedunderlying theories remain unclear. As there is no perfecttheory with which to expound the convergence and its rate, the skill of each algorithmand its modification cannot be estimatedaccurately and thereforeit is often based on many numerical tests. Nevertheless, because of the high computational cost of MEDMs and the relatively new PO research, researchon the skill of different algorithms in MEDMs is rare and their optimization effciency remains unclear.

3 SKILLASSESSMENTSTRATEGIES FOR MODEL VALIDATION 3.1 Different methods to estimate the goodness-of- fit ofoptimized models

Model validation generally analyzes and estimates the goodness-of-fit of optimized modelsto validate the effectiveness of model calibration, and it is commonly based on the skill assessment of the optimized models.However, thereis no specific method with which to validate optimized models. Commonly, modelersattempt to determinewell the model fits the observational data from a number of different aspects.The most commonmethod is visual comparison by graphical analysis (Gregg et al., 2009; Stow et al., 2009) . By setting different definitions for eachcoordinate, the goodness-of-fit can be analyzed from various respects. For example, time series plots of simulations and observations appear the common standard, as shown by Arhonditsis and Brett (2004) . Bivariate plots of simulations versus observations and thegraphical examination of residuals or misfits are also employed, as has been illustrated by Stow et al. (2009) . Additionally, comparison of spatial maps can also be adopted to test the spatial goodness-of-fit, as shown by Rose et al. (2009) .Graphical analysis is a type of qualitative description of the goodness-of-fit, while quantitative methods refer to statistical analyses and they are used mainly for the calculation of metrics such as the correlation coeffcient and root mean square error (Gregg et al., 2009; Stow et al., 2009) .

Generally, considering the computational costs, available data, modelling purpose, and others, only several state variables of interest are calibrated. However, because of interactions to differentdegrees among the state variables, only optimizing those variables of interest will affect others and result in calibration bias (Arhonditsis and Brett, 2004) .The response of unassimilated variables to PO should be analyzed to validate unassimilated variables (Gregg et al., 2009) . Comparison between the simulations and observations of all variables simultaneously considers their interactions, contributing to the analysis of the global goodness-of-fit. In contrastto single-variable methods, multivariate comparison methods such as multivariate statistical analysis and the calculation of cost function by integrating misfits of multiplevariables (Stow et al., 2009) , can better reveal the calibration effectiveness, especiallywhen multi-objective optimization is used.

Because of the comparatively low attention paid to model validation, methods for model validation tend to be mainly based on visualcomparison rather than quantitative validation. Arhonditsis and Brett (2004) reviewed 153 aquatic biogeochemical models published from 1990-2002. They found that 30.1% of the studiesreported goodness-of-fit measures, while others oftenonly presented a basic visualcomparison, usually a time series plot.A following analysis reported no relationship betweenthe level of skill assessment presentedor the accuracy of the model and the subsequent citation rate of the published paper (Arhonditsis et al., 2006) . Additionally, in the field of oceanography, Stow et al. (2009) reviewed 142 papers on marinebiological/physical models published from 2000-2007. They found that 31.7% of the papers provided quantification of the goodness-of-fit. The greater demandsof managers and decision makers mean that model validation is becoming increasingly necessary and important.

3.2 Different validation data used for skill assessment of optimized models

In addition to validation methods, different observational data used in the validation (which are called validation data) can also show different aspects of model skill. Accordingto whether the validation data are within the spatiotemporal scale covered by the observational data used in the calibration (which is called the calibration scale) , replicative or predictive performance can be assessedto describe different levels of model skill (Arhonditsis and Brett, 2004) . Replicative performance means the abilityof an optimized model to replicatethe observations within the calibration scale, depending on the skill of the optimization algorithmand the observational data used in the calibration (which are called the calibration data) . The data used in replicative calibration can be a subset of the same data set as the calibration data, or another data set within the calibration scale. The predictive performance refers to the forecasting ability of the optimizedmodel outside the calibration scale, depending on the skill of both the optimization algorithm and the model.Herein, validation data are independent of the calibration data. When the study area or physical environment of these independent data sets is different from where the optimized model is built, the performance to be validatedis also called model transferability or portability, which shows whether the structure of the model can be extrapolated to a different ecosystem (Arhonditsis and Brett, 2004; Friedrichs et al., 2007) . A model with good predictive performance that always demonstrates high reliability and robustness is desirable to support scientific research and to be the basic tool of marineenvironment management (Gregg et al., 2009) .

Replicative and predictive performances reflect two differentlevels of model skill. Models with good predictive performance are generally well replicative, but not vice versa.Overfitting and equifinalityare two typical examples.Overfitting means that misfits are arbitrarily small, such that the model is well replicative. However, this is commonlyachieved by setting free parameters, relativeto the control parameters, with values that are unrealistically too large or small that would certainlydecrease the predictive performance. Equifinality refers to model precision that remains almostthe same while the parameters are very different. This is a relatively common phenomenon in the field of watershed hydrology modelling, which also appearsin the PO of MEDMs (Ward et al., 2010) .The phenomena of overfitting and equifinalityboth make modelswell replicative but poorlypredictive, which can be diagnosed by predictive validation. The most direct cause of these phenomena is the lack of a perfect theory aboutthe termination criterion of optimization algorithms, especially stochastic global optimization algorithms. However, the primary cause is related to the variousand complex uncertainties of MEDMs and therefore, predictive validation is an important aspect of the skill assessment of MEDMs.

4 DISCUSSION 4.1 Global treatmentsrequired by the complexity of MEDMs

The marine ecosystemis complex and affected by global change and anthropogenic perturbations. Hence, MEDMs generallyconsider physical, chemical, and biological processes simultaneously, often makingnonlinear correlations of different degrees between both the parameters and the state variables. Thus, integrated globaltreatments of PO and skillassessment are requiredto recognize and to improve model skill.

Global treatments of model parameters include two aspects.First, to reduce the dimensions of the solution space of the optimization problem, parameters with large impacts on the model results but without strong inter-correlations should be selected as the set of control parameters (Chu et al., 2007; Cossarini and Solidoro, 2008; Gibson and Spitz, 2011; Zheng et al., 2012; Morris et al., 2014) . In other words, selectionof the control parameters should consider the impact of the parameters on the model results and their inter-correlationssimultaneously. The method combining LS A with the statistical approach is feasible, although its effectiveness needs to be validated (Solidoro et al., 2003; Shi et al., 2014) . Furthermore, GSA simultaneously tests the impact of multiple parameters and their interactions on the modelresults to obtain global parameter sensitivities, and it should be considered the basic method with which to select the control parameters (Cossarini and Solidoro, 2008; Gibson and Spitz, 2011; Zheng et al., 2012; Morris et al., 2014) .

The complexity of MEDMs makesthe related objective function a multimodalfunction and the optimization problemmulti-extremum. It requiresthe optimizationalgorithm to have the capability of globally searching the solution space, which has not yet been realized effectively in complex systems.The popular approach is to escape local convergence by stochastic searching and thus, to obtain an approximate or quasi-global optimum. Currently, the PO of MEDMs commonlydepends on local algorithms such as the variational adjointand Kalman filter (Jones, 1965; Evensen, 1994; Zhao et al., 2005; Fan and Lv, 2009; Gregg et al., 2009; Ward et al., 2010) . Although model precisionhas improved greatly, predictive performance might be unsatisfactory and thus, applications may be limited.The stochastic global optimization algorithm, which is generally based on the imitationof natural phenomenaand social behaviors, is being explored in the PO of MEDMs (Matear, 1995; Hurtt and Armstrong, 1999;Kavanagh et al., 2004;Oschlies and Schartau, 2005; Weber et al., 2005; Ward et al., 2010; Miró et al., 2012) . However, because it is a recent proposal, its advantages in global optimization have yet to be widely recognized. Furthermore, additional studies on the comparisons and modificationsof these stochastic globaloptimization algorithms in thefield of MEDMs are encouraged, becausethe skill assessment of algorithms is dependent on the number of numerical experiments.

In relationto the state variables of MEDMs, it has been proventhat efforts to increasemodel complexity do not necessarily lead to a systematic improvement in modelskill (Matear, 1995; Hurtt and Armstrong, 1999; Arhonditsis and Brett, 2004;Flynn, 2005; Friedrichs et al., 2007, 2009; Xiao and Friedrichs, 2014a) . Thisis mainly becauseavailable observational data are generally insuffcient to support models with relatively high complexity. Thus, the number of state variables of most MEDMs is relatively small (generally<10) and much smaller than the number of parameters (usuallydozens or hundreds) . Compared with comprehensive treatmentson multiple parameters, it is also relatively effective to consider only a single state variable at a time. However, as the number of state variablesincreases, multi-objective optimization and multivariate comparison become essential.

4.2 Effcient use of available data to direct model skill improvement

Observational data are the important part of modelling. It is more effcient to undertakemodel calibrationor validation at multiple sites simultaneously, whereas in most cases, in situ data are unavailable on the spatiotemporal scale of interest (Xiao and Friedrichs, 2014b) .This is one reason why research on the skill of MEDMs develops slowly (Arhonditsis and Brett, 2004; Gregg et al., 2009) . Fortunately, satellite-derived data of the marine ecosystem have brought significant progress during the last two decades (Rousseaux et al., 2013; Xiao and Friedrichs, 2014b) . Despite the limitation of only providing information in the first opticaldepth, satellite-derived data provide comprehensive synoptic coverage over large regions of the ocean, the advantages of which in terms of skill improvement have been illustrated by Xiao and Friedrichs (2014b) . Furthermore, the problem of how to exploit fully the limited observational data, to revealand improve model skill, remains an importantissue.

In additionto the necessary data used for the physical forcing, initial conditions, and boundary conditions of each state variable, model calibration requires at least one data set. When thereis only one data set, if it is used completely in calibration to pursue the maximum improvement in model precision, there are no available data for validation to quantify the skill of the optimized model.Therefore, a subset shouldbe retained for replicative validation. In such a case, the sizes of the two subsetsshould be allocated rationally to guarantee the effectiveness of both calibration and validation. Gregg et al. (2009) suggested to first apply the entire data set to the calibration to achieve a set of hopefully optimal parameters. Then, the data should be divided into two equal subsetsfor calibration and validation purposes, to observesimultaneously the deterioration of the calibration relativeto the improvement in the validation.

When there are two sets of independent data, one is usually used for calibration and the other for validation. Thus for model validation, modelers should measurevarious metrics and combine graphical analysis with quantitative methodsto reveal comprehensively the skill of the optimizedmodel. In addition, using the two data sets flexiblymeans diverse strategies for calibration can be used to optimize and assess model skill. For example, the two data sets could both be appliedto model calibration to obtain two relativesets of optimalparameters. Then, one data set could be used to validatethe relative parameter sets of the other, which is generally called cross validation (Friedrichs et al., 2009) . Moreover, by combining the two data sets for model calibration, a common set of optimalparameters for different study areas or ecosystems can be obtained. Additionally, learning about the climaticconditions and physicalenvironment under which the available data are observed and analyzing the observational data mathematically are necessary for modelers to understand the real worldthat is to be modeled, which can direct modelersto design the model structure scientifically.

Generally, observational data are availableto modelers to assistin the improvement of model skill; however, in cases where there is a lack of available observational data, approaches using simulated data are adopted.A typical methodology is the parameter SA, which analyzes the degree of impact of parameters using simulation data with different values of the parameters (Fennel et al., 2001;Solidoro et al., 2003; Shi et al., 2012, 2014) . Additionally, comparison of simulations resulting from an optimized model and free-run model is also useful when combined with observations to test the ability of differentoptimization algorithms to improvemodel precision. Regarding DA, a typical application is to use simulation data for calibration, which is calledan identical twin experiment (Friedrichs, 2002) . It can be used to compare and test the skill of different optimization algorithms withoutobservational data, and it can help in the examination of the effect of calibration data on optimization effciency; thus, guiding the design of the samplingof the observational data (Friedrichs, 2002) .

4.3 Model uncertainty analysisfor impact factors on model skill

Current techniques used for the improvement in model skill have been able to achieve the requirement of replicating real natural systems.Nevertheless, managers and decisionmakers pay considerable attention to the reliability and robustness of models during the process of model application, which places high requirements on predictive validation (Gregg et al., 2009; Fulton, 2010; Ward et al., 2010) . MEDMs are affected by various factorswhose errors contribute to model uncertainties, which results in prediction uncertainty. By analyzingthe impact of different factors on the model results, model uncertainty analysis (UA) can supply model users with a quantitative assessment of model confidence, which can assistin model applications and is an important aspect of model skill assessment.

According to the impactfactors on modelskill, model uncertainties are derived from four aspects: model parameters, structure, solution, and inputs (Gregg et al., 2009) . Model parameter uncertainty is derived from uncertainties of the PO algorithm and calibration data, as well as the inherent variations of parameters in time and space. Model structure uncertainty is introduced both from insuffcient knowledge of the actual problems and from treatments of conceptualization and simplification on natural processes. Model solution uncertainty means the approximate errors of the numerical algorithm. Model inputs uncertainty is derived from both the observational errors and the randomness of natural processes. Existingstudies have shownthat uncertaintiesof MEDMs are caused mainly by observational errors and model parameters. Therefore, UA of MEDMs focuseson the uncertainties of the inputs and parameters, as well as on the prediction uncertainty caused by the various factors in combination.

Currently, researchof UA in MEDMs overlaps with model calibration and validation (Wallhead et al., 2009) . For example, using the Sobol method, the contribution ratesof each parameter and their interactions to model uncertainties are used to describe model parameter uncertainty (Gibson and Spitz, 2011; Morris et al., 2014) .In addition, according to the idea of SA, model prediction uncertainty caused by the inputs is analyzedby applying disturbances with different degreesto the input data (Friedrichs et al., 2009) . The Bayesian approach, which is distinguished fromthe traditional frequency approach in its treatmentof the model inputs as random variables, is considered a good choice for model UA (Ward et al., 2010) . It includes the earliest methodof GeneralizedLikelihood Uncertainty Estimation (Beven and Binley, 1992) and the widely used approach of the MarkovChain Monte Carlo (Kuezera and Parent, 1998) . It has been commonly used in the fields of watershedhydrology, water quality modelling, and oceanography, but mainlywith simple models. For example, Stow and Scavia (2009) implemented a Bayesian hierarchical version of a simple process-based model to researchhypoxia in Chesapeake Bay, and they quantified the prediction uncertainty caused by parameters and inputs using an ensemble estimat ion. Combininga priori and a posteriori information will help direct an integrated analysis of model uncertainties to reduce prediction uncertainty (Ward et al., 2010) . Therefore, the application of the Bayesianapproach in MEDMs should be explored further.

5 CONCLUSION

Based on a technical framework (shown in Fig. 1) , the methods andstrategies used to improve and assess the skill of MEDMs were discussed. The objective was to inspire further ideas regarding the skill improvement of MEDMs and to establish a standard technical framework for its continuedscientific development. Particularly, global treatments on SA and PO were emphasized. A feasible strategyfor GSA is first to excludethose parameters that have the least impact on the model resultsusing the Morrismethod, and then to quantify the contribution rates of the uncertainties of those parameterswith significant impacts using the Sobol method. Moreover, stochastic global optimization methods (shownin Table 1) should be considered seriouslyduring PO, of which the micro-genetic algorithmand simulated annealing methods are recommended. It should be noted that the necessity for in-depth researchof the model skillof MEDMs depends on the modelling purpose, rather than modelers’interest and ability. However, of course the efforts made in improving model skill, described in this paper, are not intendedto block further researchprogress and paper publication. The actual research should consider equally scientific advancementand management practice.