Adaptive Time Step Algorithms for the Simulation of marine Ecosystem Models using the Transport Matrix Method Implementation Metos3D

The reduction of the computational effort is desirable for the simulation of marine ecosystem models. Using a marine ecosystem model, the assessment and the validation of annual periodic solutions (i.e., steady annual cycles) against observational data are crucial to identify biogeochemical processes, which, for example, influence the global carbon cycle. For marine ecosystem models, the transport matrix method (TMM) already lowers the runtime of the simulation significantly and enables the application of larger time steps straightforwardly. However, the selection of an appropriate time step is a challenging compromise between accuracy and shortening the runtime. Using an automatic time step adjustment during the computation of a steady annual cycle with the TMM, we present in this paper different algorithms applying either an adaptive step size control or decreasing time steps in order to use the time step always as large as possible without any manual selection. For these methods and a variety of marine ecosystem models of different complexity, the accuracy of the computed steady annual cycle achieved the same accuracy as solutions obtained with a fixed time step. Depending on the complexity of the marine ecosystem model, the application of the methods shortened the runtime significantly. Due to the certain overhead of the adaptive method, the computational effort may be higher in special cases using the adaptive step size control. The presented methods represent computational efficient methods for the simulation of marine ecosystem models using the TMM but without any manual selection of the time step.


Model Description
A marine ecosystem model represents the interaction between the ocean circulation, the ocean biota and marine biogeochemical cycles. This involves modeling the marine ecosystem by a given number of ecosystem species (or biogeochemical tracers), which are substances in the ocean water and subject to chemical or biochemical reactions. Due to the full coupling of the ocean circulation with the tracers, in which the circulation influences the tracer concentrations and, vice versa, the tracer concentrations affect the circulation, the simulation of a fully coupled model (online model) is computationally expensive. Restricting this coupling to the influence of the ocean circulation on the tracer concentrations (i.e., an offline model using passive tracers that do not affect the ocean circulation) reduces the computational effort. In particular, this one-way coupling enables the application of a pre-computed ocean circulation for the simulation.

Model Equations for marine Ecosystems
A system of partial differential equations describes the marine ecosystem model. The number of modeled tracers defines the complexity of the marine ecosystem model and, hence, the size of the system of differential equations. For n y ∈ N tracers on a spatial domain Ω ⊂ R 3 (i.e., the whole global ocean) and a time interval [0, 1] (i.e., one model year), we consider marine ecosystem models using an offline model. With the function y i : Ω × [0, 1] → R, i ∈ {1, . . . , n y }, of the tracer concentrations for the single tracer y i and the vector y := (y i ) ny i=1 of all tracers, the system of parabolic partial differential equations i ∈ {1, . . . , n y }, represents the tracer transport of a marine ecosystem model. In Eq. (1), D denotes the diffusion, A the advection operator, and q represents the biogeochemical model describing the interaction between the different tracers. The homogeneous Neumann boundary conditions (2) including the normal derivative model that there are no fluxes on the boundary. The advection and diffusion, coming from the ocean circulation, determines the tracer transport in the ocean. For a given velocity field v : Ω × [0, 1] → R 3 , the linear operator A : Ω × [0, 1] → R models the advection as A(x, t)y i (x, t) := div (v(x, t)y i (x, t)) , x ∈ Ω, t ∈ [0, 1], for i ∈ {1, . . . , n y }. The diffusion D : Ω × [0, 1] → R used as model of the turbulent effects of the ocean circulation neglects the molecular diffusion of the tracers themselves since this is known to be much smaller in relation to the diffusion caused by turbulence. As a result of the quite different scales in horizontal and vertical direction, the diffusion operator is split into a sum D = D h + D v of a horizontal and vertical part in order to treat the vertical part in the time integration implicitly. The diffusion is split into horizontal and vertical part as for i ∈ {1, . . . , n y }. Here, div h and ∇ h denote the horizontal divergence and gradient, respectively, and z the vertical coordinate. The diffusion coefficient fields κ h , κ v : Ω × [0, 1] → R are the same for all tracers due to the fact that the molecular diffusion is generally assumed to be smaller than the diffusion induced by the turbulence of the ocean circulation.
The biogeochemical model contains the biogeochemical processes within the marine ecosystem. The interplay of the biogeochemical model with the effects of the ocean circulation (i.e., the whole system (1) to (6)), on the other hand, is called the marine ecosystem model. For the tracer y i , i ∈ {1, . . . , n y }, the, in general, nonlinear function q i : Ω × [0, 1] → R, (x, t) → q i (x, t, y, u) describes the biogeochemical processes for this tracer. In particular, this nonlinear function q i includes firstly the influence of the variability of the solar radiation at space and time, secondly the coupling of this tracer to the other species and thirdly n u ∈ N model parameters u ∈ R nu controlling, for example, growth, loss and mortality rates or sinking speed of this tracer. Altogether, the biogeochemical model q = (q i ) ny i=1 summarizes the biogeochemical processes of all tracers. An annual periodic solution of the marine ecosystem (i.e., a steady annual cycle) satisfies in addition to (1) and (2) y i (x, 0) = y i (x, 1), x ∈ Ω, (6) for i ∈ {1, . . . , n y }. For this purpose, we assume that the operators D, A and the functions q i are also annually periodic in time.

Semi-discrete Setting
For the adaptive step size control, we used a semi-discrete setting where the computational domain is already discretized in space, but time is kept continuous. Then Eq. (1), (2) read ∂ y i ∂t (t) + (D(t) + A(t)) y i (t) = q i (t, y(t), u), t ∈ [0, 1], with initial value y i (0) = y 0 i .

Biogeochemical Models
The biogeochemical models, used in this paper, differed in the number of ecosystem species. [15] introduced a hierarchy of five biogeochemical models with an increasing complexity starting with a simple model including only one tracer up to a model with five tracers. In addition to this hierarchy, we applied the biogeochemical model of [7]. In the following, we briefly describe the biogeochemical models and refer the reader to [15,27,7] for a detailed description of the modeled processes and model equations. Table 1 summarizes the model parameters of the different biogeochemical models. Especially, the biological production depends on the available light. The light intensity decreases with the depth wherefore the ocean is divided into two layers, a euphotic (sun lit) layer of about 100m and an aphotic zone below. Depending on the insolation based on the astronomical formula of [22], the light limitation function I : Ω × [0, 1] → R ≥0 models the available light taking the ice cover and the exponential attenuation of water into account. Since the main part of the biological production occurs in the euphotic layer, particulate matter sinks from the euphotic layer to depth and remineralizes there according to the empirical power-law relationship [19].
The N model is the simplest biogeochemical model of the hierarchy and represents only phosphate (PO 4 ) as inorganic nutrients, cf. [1,15]. The available nutrients and light restrict the phytoplankton production (or biological uptake). The phytoplankton production depends on the maximum production rate µ P and applies an implicitly prescribed concentration of phytoplankton y * P = 0.0028 mmol P m −3 . Altogether, n u = 5 model parameters listed in Table 2 control the biogeochemical processes of the nutrient tracer y = ( y N ).
Implicit representation of sinking speed 1 a D Increase of sinking speed with depth The N-DOP model includes dissolved organic phosphorus (DOP) in addition to nutrients (N), i.e., y = ( y N , y DOP ), cf. [2,23,15]. This model computes the phytoplankton production also with (8) and introduces n u = 7 model parameters ( Table 2).
The NP-DOP model contains phytoplankton (P) in addition to nutrients (N) and dissolved organic phosphorus (DOP), i.e., y = ( y N , y P , y DOP ), cf. [15]. As a result of the explicit treatment of phytoplankton, the NP-DOP model computes the phytoplankton production again with (8) but using the explicit phytoplankton concentration y P instead of y * P . A quadratic loss term for phytoplankton, moreover, models the zooplankton grazing with the implicitly prescribed zooplankton concentration y * Z = 0.01 mmol P m −3 . The n u = 13 model parameters listed in Table 2 control the biogeochemical processes of the NP-DOP model.
The NPZ-DOP model consists of four tracers, nutrients (N), phytoplankton (P), zooplankton (Z) and dissolved organic phosphorus (DOP), i.e., y = ( y N , y P , y Z , y DOP ), cf. [15]. While the phytoplankton production (8) is the same as for the NP-DOP model, the zooplankton grazing (9) explicitly contains the zooplankton concentration y Z instead of the implicitly prescribed concentration y * Z . The model introduces n u = 16 model parameters summarized in Table 2.
The NPZD-DOP model, the most complex model of the hierarchy, finally, contains detritus (D) in addition to nutrients (N), phytoplankton (P), zooplankton (Z) and dissolved organic phosphorus (DOP), i.e., y = ( y N , y P , y Z , y D , y DOP ), cf. [34,15]. Both the phytoplankton production (8) and the zooplankton grazing (9) are identical to the NPZ-DOP model. Table 2 lists the n u = 18 model parameters of the NPZD-DOP model.
The MITgcm-PO4-DOP model contains phosphate (PO 4 ) and dissolved organic phosphorus (DOP). This  [7] resembles the N-DOP model. It also has n u = 7 model parameters which we identified with the models parameters of the N-DOP model ( Table 2).

Transport Matrix Method
The transport matrix method (TMM) reduces the simulation of the tracer transport of an offline model to matrixvector multiplications. [13] applied a linear matrix equation instead of directly implementing a discretization scheme for the advection-diffusion equation (1) because the application of the operators A and D for the advection and diffusion on a spatially discretized tracer vector is linear, see also [11]. Consequently, the TMM approximates the ocean circulation by matrices which include the influence of all parameterized processes represented in the underlying ocean circulation model on the transport. Each time step of the simulation with the TMM requires only two matrix-vector multiplications and an evaluation of the biogeochemical model. For the discretization of the advection-diffusion equation, let (x k ) nx k=1 a spatial discretization with n x ∈ N grid points of the domain Ω (i.e., the ocean) and the time steps t 0 , . . . , t nt ∈ [0, 1], n t ∈ N, specified by t j := j∆t, j = 0, . . . , n t , ∆t : an equidistant grid of the time interval [0, 1] (i.e., one model year). For a time instant t j , j ∈ {0, . . . , n t − 1}, the vector y ji ≈ (y i (t j , x k )) nx k=1 ∈ R nx is a spatial discretization of the tracer y i , i ∈ {1, . . . , n y }, and q ji ≈ (q i (x k , t j , y j , u)) nx k=1 ∈ R nx the spatially discretized biogeochemical term q i for the tracer y i . Besides, y j := ( y ji ) ny i=1 ∈ R nynx and q j := ( q ji ) ny i=1 ∈ R nynx summarize the tracer discretization as well as the biogeochemical terms for all tracers using a reasonable concatenation. Discretizing the advection and horizontal diffusion explicitly and the vertical diffusion implicitly, the application of a semi-discrete Euler scheme for (1) yields a time-stepping y j+1 = I + ∆tA j + ∆tD h j y j + ∆tD v j y j+1 + ∆t q j ( y j , u) , j = 0, . . . , n t − 1, with the identity matrix I ∈ R nx×nx and the spatially discretized counterparts A j , D h j and D v j of the operators A, D h and D v at time instant t j , j ∈ {0, . . . , n t − 1}. The matrices D v j are block-diagonal since they involve only the vertical part of the diffusion. Thus, each water column is separated from the others. We denote the explicit and implicit transport matrices by for each time instant t j , j ∈ {0, . . . , n t − 1}. The implicit matrices T imp j , j ∈ {0, . . . , n t − 1}, are block-diagonal again since the inversion of a matrix keeps this structure. Finally, a time step of the marine ecosystem model using the TMM is given by Due to the grid-point based ocean circulation model, both the explicit and the implicit transport matrices are sparse. For the implicit ones, this is given by their block-diagonal structure.
In practical computations, the TMM computed and stored monthly averaged matrices and interpolated those linearly for any time instant t j , j = 0, . . . , n t − 1. In this paper, we applied transport matrices computed with the MIT ocean model [18] used a global configuration with a latitudinal and longitudinal resolution of 2.8125 • and 15 vertical layers.

Computation of steady annual Cycles
For a marine ecosystem model, the steady annual cycle is a fixed-point of the spin-up. Applying the above iteration (14) over one model year, the steady annual cycle (i.e., an annual periodic solution) in a fully discrete setting fulfills The steady annual cycle is a fixed-point of the nonlinear mapping with ϕ j defined in (14), describing the time integration of (14) over one model year. A classical fixed-point iteration takes the form using an arbitrary start vector y 0 ∈ R nynx and model parameters u ∈ R nu . Interpreting the fixed-point iteration as pseudo-time stepping or spin-up, the vector y ∈ R nynx contains the tracer concentrations at the first time instant of the model year ∈ N. We tested the numerical convergence of the spin-up (i.e,. of the iteration (17)) with the difference between two consecutive iterates determined by ε := y − y −1 (18) for iteration (model year) ∈ N 0 . For this purpose, we quantified this difference with various norms. Namely, we defined a weighted Euclidean norm with weights w k ∈ R >0 for k ∈ {1, . . . , n x } and z ∈ R nynx indexed as . If all weights are equal to 1 (i.e., w k = 1 for k = 1, . . . , n x ), the norm · 2,w corresponds to the Euclidean norm · 2 . We denoted by · 2,V the discretized counterpart of the L 2 (Ω) ny norm using the weights w k = |V k |, k = 1, . . . , n x , with the box volume |V k | corresponding to the grid point x k . In order to consider not only the difference at first time instant but for the whole trajectory over one model year, we, moreover, defined the weighted Euclidean norm · 2,w,T by . Analogous to the weighted Euclidean norm · 2,w , we denoted by · 2,T the Euclidean norm and by · 2,V,T the weighted Euclidean norm · 2,w,T using weights w k = |V k | for k ∈ {1, . . . , n x }. In addition to the norms including all grid points of the discretization (x k ) nx k=1 , we restricted norms to horizontal layers of the ocean discretization, such as to the upper layer describing the ocean surface. We identified with · | L for L ⊂ {1, . . . , 15} the restriction of norm · to the layers selected in the set L. For example, · | L with L := {1} restricted the norm to the grid points of the ocean surface.

Temporal Coarsening of Transport Matrices
Using simple matrix operations, transport matrices computed with a given time step can be used to generate matrices corresponding to a bigger time step. The procedure has been described in [11]. Following his approach, we used the transport matrices T exp j and T imp j , j ∈ {0, . . . , n t − 1}, to generate transport matrices corresponding to a time step with coarsening factor m ∈ N. Consequently, these matrices represent larger time steps than the ones in the underlying ocean circulation model which were used to generate the transport matrices T exp j and T imp j . The explicit transport matrix T exp j,m is the exact representation of the larger time step, i.e., The implicit transport matrix T imp j,m is an approximation with a loss of accuracy of order ∆t 2 since (using the binomial theorem) cf. [28,11]. These transport matrices are still sparse [11]. The time step of the ocean circulation model used for the computation of the transport matrices T exp j and T imp j , j ∈ {0, . . . , n t − 1}, corresponds to 3 h [18]. Assuming 360 days per model year, the number of time steps per model year is n t = 2880. We, hereinafter, denoted this time step with 1 ∆t using ∆t = 1 2880 . More specifically, time step 2 ∆t corresponds to the coarsening factor 2, i.e., a doubling of the effective time step, with n t = 1440. In order to identify the time step used to run the spin-up (17), we denoted by the time-integration of one model year using the time step m∆t with the factor m ∈ N for y ∈ R nynx and u ∈ R nu . More importantly, the transport matrices T exp j,m and T imp j,m , j ∈ {0, . . . , n t − 1}, enter (14) for the application of the time step m∆t instead of the transport matrices T exp j and T imp j . Any choice m ∈ N is possible in the above computations. For every value of m to be used, twelve pairs of explicit and implicit transport matrices, however, have to be supplied. Since we wanted to check a rather wide variability of the time steps, we have chosen the possible coarsening factors m ∈ M := {1, 2, 4, 8, 16, 32, 64} (28) in this study. The upper limit of m = 64 was motivated by our observation that the spin-up for most of the models did not converge anymore for larger time steps.

Negative Tracer Concentrations
A time step that is too large can cause negative tracer concentrations. Given a non-negative tracer distribution y j in all grid points at a time instant t j , j ∈ {0, . . . , n t − 1}, the multiplication with the explicit transport matrix T exp j in (14) will always result in a non-negative distribution. In contrast, the source term q j ( y j , u) may give negative values for some tracers at some points even though the input variable y j is non-negative. On the other hand, negative tracer concentrations do not make sense. Moreover and without going into mathematical details, it can be shown that at least in the semi-discrete, time-continuous setting, a non-negative initial value always leads to a non-negative tracer distribution for all time, at least for the biogeochemical models used here. Thus, negative values may be obtained only due to a time step that is too big. Following this idea, one might consider to guarantee non-negativity throughout the simulation, in order to keep the solution consistent with biogeochemistry and mathematics. In many models, non-negativity is enforced by a simple setting of negative values to zero before the evaluation of the source term. Such a setting obviously increases the total mass in the ecosystem. Hence, a frequent occurrence of big negative values and their correction to zero will result in a changed steady solution obtained in the spin-up. Despite such a correction before the source term step, the sum in the bracket in (14) might be negative at some points, depending on the used time step. Thus, after applying the implicit matrix, the result of the time step y j+1 might contain negative values as well.
Small negative values in some points do not necessarily change the convergence of a spin-up to a steady annual cycle. Therefore, a main criterion for algorithms that use bigger time steps (as the ones presented here) is if the spin-up converges, and if it converges to the same solution as for the standard setting. We also investigated an algorithm that enforces non-negative tracer concentrations and how this effects the solution and the cost saving.

Step Size Control Algorithms
Methods for automatic step size control automatically adapt the time step during the calculation of a transient computation. They may, therefore, also be used for the steady annual cycle computation via a spin-up. A step size control method estimates the local discretization error by computing two approximations, either with different time steps (in the method used here) or by using two different time integration methods. The approach used here is based on the Richardson extrapolation [31,30]. For the estimate of the local discretization error, the step size control computes two approximations with two different time step sizes in every step. Depending on the error estimate and a desired accuracy, the step size control accepts or rejects the approximation calculated with the smaller time step (the approximation computed with the larger time step always serves the error estimation only). Then, it adapts the step size to a value that, using the estimation of the error, would result in the desired accuracy. Thus, in this step, an increase or a decrease of the step size is possible. Using the adapted step size, the method, afterwards, either starts the calculation of the next step or reruns the calculation of the current step. In summary, the step size control finally uses always the largest time step that keeps the error below the given tolerance.
Step size control methods are based on the existence of an asymptotic expansion of the discretization error [37,Sect. 7.2.3]. To obtain this result, the unknown solution is assumed to be three times continuously differentiable in time. The question whether this regularity is satisfied for the models investigated in this work is not discussed here. We just note that all models have continuous right-hand sides (i.e., source-minussink terms). Furthermore, it is quite usual to apply a step size control even though some assumptions of the underlying mathematical theory may not be given or shown for an actual application. Since we compared our results obtained with the step size control to those with constant time steps, an assessment of the method is possible.
Algorithm 1 depicts the step size control for the calculation of a steady annual cycle. It is based on [6, Algorithm 5.2]. Since we are aiming at a steady annual cycle obtained in a spin-up, we do not check for the local discretization error after every single time step, but only after an adjustable number of model years. This number (n s ∈ N) is one of the input parameters of the algorithm. The loop over the n s years is realized in lines 10 to 13. Another important feature is the set of possible or desired multiplication factors m. Here, the value m = max M = 64 is inadmissible because the step size control estimates the local error using two approximations calculated with the selected time step m∆t and a larger time step. Thus, the setM := {1, 2, 4, 8, 16, 32} defines the admissible time steps of Algorithm 1.
We designed the step size control algorithm having various options. The user may select, firstly, the initial time step m init ∆t of the step size control with m init ∈M . The second option is the number n s ∈ N of model years after which the error is estimated. The step size control computes the two approximations with time steps m∆t andm∆t for the model year + n s , ∈ N 0 , and m,m ∈ M with m <m and estimates the local discretization error y +ns,m − y +ns,m . The choice of the next time step depends on this error (line 15 or 21 and following). Thirdly, it is possible to choose the norm that may have an influence on the error estimation. Another setting of the step size control, fourthly, is the tolerance τ 0 ∈ R >0 which controls the acceptance of the approximation. By default, the step size control applies the coarsening factor m init = 1 for the initial time step m init ∆t, n s = 1, the volume-weighted Euclidean norm · 2,V and the tolerance τ 0 = 1.
A variant of Algorithm 1 includes an additional avoidance of negative concentrations (see Algorithm 2). If enabled, the step size control accepts an approximation -as additional criterion to the local discretization error

Decreasing Time Steps Algorithm
The step size control algorithms described above require an increased computational effort to automatically adjust the time step. They both try to find the optimal step size in the sense that it should be as big as possible but still keeps the local discretization error below the given tolerance. Therefore, the algorithms are automatically able to both increase and decrease the time step during the spin-up. This optimality comes with the additional effort of computing always two approximations to estimate this error.
The decreasing time steps algorithm presented in this section, in contrast, exclusively decreases the time step during the spin-up. The motivation is the same, namely to reduce the computational costs. The algorithm assumes that in its beginning bigger time steps are sufficient. The procedure is shown in Algorithm 3. Starting

Algorithm 2
Step size control avoiding negative concentrations. Shown are only differences to Algorithm 1. Input: Initial concentration: y 0 ∈ R nynx , Parameter vector: u ∈ R nu , Number of model years: T ∈ N, Coarsening factor of the initial time step: m init ∈ M , Number of model years for the error estimation: n s ∈ N, Tolerance: τ 0 ∈ R >0 Output: y T ∈ R nynx // Lines 1 to 14 identical as in Algorithm 1  15 = + n s 12: end while from an initial time step m init ∆t with a coarsening factor m init ∈ M , the algorithm uses the same time step until the norm y − y +ns , ∈ N 0 , fell below a given tolerance ε ∈ R >0 . In other words, the spin-up calculation with this time step already reached almost a steady annual cycle. If the norm fell below the tolerance, the algorithm decreases the time step. The assessment whether a significant reduction was still achieved with the current time step takes place after a fixed number of n s ∈ N model years. By default, the decreasing time steps algorithm uses the initial time step m init ∆t with the coarsening factor m init = 64, n s = 50 and the tolerance ε = 0.001.

Results
The automatic adjustment of the time step used for the spin-up computation influenced both the computational effort and the accuracy of the approximation of the steady annual cycle. In this section, we present the numerical results using the adaptive step size control (see Sect. 3) and the decreasing time steps algorithm (Sect. 4) to shorten the runtime of the computation of the steady annual cycle. We assessed the accuracy and cost saving of the calculated approximations.

Experimental Setup
For each part of the spin-up calculation (i.e., n s ∈ N model years), we used the marine ecosystem toolkit for optimization and simulation in 3D, Metos3D, [27,28]. Overall, we computed all spin-ups over 10000 model years. We started the spin-up always with a constant global mean tracer concentration of 2.17 mmol P m −3 for PO 4 and, if present, 0.0001 mmol P m −3 for the tracers DOP, P, Z and D.
Using different parameter vectors, we calculated the steady annual cycles for the various biogeochemical models. We used, on the one hand, the reference parameter vectors listed in Table 2 and, on the other hand, 100 parameter vectors generated by a Latin hypercube sample within the bounds defined in Table 2 for each biogeochemical model, cf. [20]. We created these parameter vectors by the lhs routine of [17].
For the assessment of the calculated approximations of the steady annual cycle, we compared them with reference solutions. These were chosen as the result obtained by a spin-up using Metos3D with constant time step m ∆t, m ∈ M , also over 10000 model years. These solutions are denoted by y 10000,m and described in more detail in [25]. In most cases, we considered the case m = 1 only. The reference solutions were also used to measure, in particular, the accuracy of an approximation x ∈ R nynx obtained by one of the algorithms (i.e., the step size control or the decreasing time steps algorithm) by the relative difference We call this quantity (29) the (relative) error of the respective result x. Furthermore, we quantified the saving of the computational costs by

Step Size Control Algorithms
In this section, we present the results using the step size control to compute a steady annual cycle (Algorithms 1 and 2). The step size control applies different time steps to adapt the step size using an estimation of the local discretization error. Starting with a default setting for Algorithm 1 (Sect. 5.2.1), we also analyzed the behavior using the avoidance of negative tracer concentrations using Algorithm 2 (Sect. 5.2.2) and, briefly, other configuration settings (Sect. 5.2.3) defined in Sect. 3.

Algorithm 1 with default Step Size Control Setting
For the reference parameter vectors of Table 2, the step size control algorithm with default setting computed approximations of the steady annual cycle that were almost identical to the ones obtained with a constant step size. For the N, N-DOP and MITgcm-PO4-DOP models, the step size control utilized the largest possible time step by increasing the time step in the first steps to 32 ∆t. Accordingly, the algorithm lowered the computational costs by 95%. In contrast, the step size control used always time step 1 ∆t (i.e., no increase of the time step occurred) for the other three biogeochemical models NP-DOP, NPZ-DOP and NPZD-DOP. Therefore, the computational effort for these three models was greater than for the solution y 10000,1 . The value computed with formula (30) was -0.5. Obviously, the obtained solution was the same as the one using constantly 1 ∆t. As an  Table 2) and the N model using the step size control with default setting (Algorithm 1) and with avoidance of negative concentrations (Algorithm 2). Shown are the convergence of the spin-up (left), i.e., the norms of the difference (18) between consecutive iterations, and the relative error (29) (right). On the left, the curves for the reference run and the one with step size control avoiding negative concentrations overlap. example, Fig. 1 indicates the similar convergence behavior towards a steady annual cycle using the step size control and the constant step size 1 ∆t, both for the N model. The convergence behavior using the step size control concurred with the one of the spin-up with constant time step 32 ∆t, i.e. y 10000,32 . In this case, the step size control doubled the time step five times up to the maximal time step 32 ∆t in the first five model years. Afterwards, the time step remained unchanged for the entire spin-up, cf. [25, Figs. 1 and 4]. Consequently, the accuracy of this approximation resembled that of y 10000,32 (Fig. 1, right). The step size control resulted in a cost saving of 95% compared to the solution y 10000,1 . Indeed, the cost saving was even 97% using directly the spin-up with time step 32 ∆t because the step size control necessitated two approximations calculated with different time steps to estimate the local error for every model year. Looking at the reduction, it would be possible to terminate the step size control even after a much lower number of model years to lower further the computational effort (Fig. 1). For the N-DOP and MITgcm-PO4-DOP model, the results were similar.
Also for the 100 parameter vectors of the Latin hypercube sample, Algorithm 1 yielded solutions that showed the same accuracy as solutions obtained with a fixed step size. For the N, N-DOP and MITgcm-PO4-DOP models, the results were consistent with those for the reference parameter vectors above, i.e., the accuracy of the approximations coincided with that for the spin-up calculated with constant time step 32 ∆t (Figs. 2 and 3). The algorithm increased the time step up to the maximum of 32 ∆t and, afterwards, used only this time step. For a few parameter vectors, the step size control, however, did not increase the time step directly, but computed several model years before increasing the time step. For these three models, the step size control resulted in a tremendous cost saving (Fig. 4). For the NP-DOP, NPZ-DOP and NPZD-DOP models and more than 85 parameter vectors, the algorithm did not increase the time step. Thus, the approximations corresponded to the reference solution y 10000,1 (Fig. 2) in these cases. For the remaining parameter vectors, at least one increase and, possibly, some decreases of the time step took place. This resulted in a larger relative error but still a reasonable accuracy of the calculated approximation. Except for a few parameter vectors where the step size control increased the time step, the step size control caused a larger computational effort than using 1 ∆t constantly because the algorithm needs two approximations to estimate the local discretization error (Fig. 4).

Relative error
Step size control Incl. avoidance Step size control Incl. avoidance  step size control avoiding negative concentrations, three instead of two approximations, hence, were computed for each model year because the step size control accepted in every step the approximation calculated with the   minimal time step 1 ∆t (independent of negative concentrations) and increased based on the error estimation the time step for the next step. In the next step, the algorithm discarded the approximation calculated with the larger time step due to negative concentrations, and reran the approximation with the minimal time step 1 ∆t again. With an adjustment of our implementation (firstly, the calculation of an approximation using the small time step to check it for negative concentrations and, secondly, the computation of the second approximation using the larger time step to estimate the local error), the computational effort could slightly be reduced to an increase to about 50%.
For the Latin hypercube sample, the same qualitative results were obtained. For the N model, the results mostly coincided with the results of the step size control using default settings (Fig. 2). Only for a few parameter vectors, the step size control reduced the time step. This was again a result of negative concentrations during the simulation. As a consequence, the accuracy improved but at the expense of the cost saving (Fig. 4). For the other biogeochemical models, the accuracy of the approximations in Fig. 2 indicates that the Algorithm 2 computed the reference solution except for some outliers because the step size control did not increased the time step during the simulation. Thus, the computational effort was larger than for the spin-up computation of the reference solution with 1 ∆t (Fig. 4).

Step size control settings
We investigated the influence of different settings for the step size control. For that purpose, we used only the reference parameter vectors listed in Table 2.
The number of model years n s ∈ N, after which the error estimation was computed, had no influence on the accuracy of the steady annual cycle approximation using the step size control (Table 3). For the N, N-DOP and MITgcm-PO4-DOP models, the step size control increased the time step as quickly as possible to the maximal time step regardless of the number of model years n s and, afterwards, applied only this maximal time step. As  a consequence, the approximations for the different model years n s were almost identical. Similar to the results above for the NP-DOP, NPZ-DOP and NPZD-DOP models, the approximation corresponded to the reference solution for the different model years n s since the step size control used always time step 1 ∆t. As can be seen in Table 4, the total cost saving depends on the parameter n s . Taking a bigger value, the algorithm needs more model years to reach the final, maximal step size. Thus, the cost saving decreases with growing values of n s . The norm used to estimate the local discretization error affected the approximation calculated with the step size control (Table 5) in some cases. Firstly, the application of the Euclidean norm · 2 instead of the norm · 2,V influenced the step size control for the N-DOP model only. For this model, the step size control, occasionally, did not increase the time step directly to the maximal time step and finished the simulation with time step 16 ∆t. This led to a smaller relative error. The use of the norms · 2,T and · 2,V,T , secondly, resulted in at least one increase of the time step for the NP-DOP, NPZ-DOP and NPZD-DOP models and, therefore, in a reduction of the computational costs. However, these norms had no influence for the N, N-DOP and MITgcm-PO4-DOP models. Lastly, the restriction of the norm · 2,V to different layers had an impact on the error estimation because the annual variability exclusively affects concentrations in the upper ocean, while the concentrations in the deeper ocean change very slowly. Besides, the box volumes in the upper ocean are small in comparison to the box volumes in the deeper ocean, but big boxes have a greater effect on the norm as small ones for a volume weighted norm. However, the local discretization error was high in the upper ocean due to the variability, and small in the deeper ocean. Thus, the step size control increased more frequently the step size using the restricted norms including deeper layers or including more layers ( Table 5).
The choice of the initial time step especially affected the three most complex biogeochemical models. Due to the directly increasing time step up to the maximum, no differences were visible for the N, N-DOP and MITgcm-PO4-DOP models ( Table 6). In contrast, the relative errors in Table 6 indicate the use of larger time steps for at least one model year for the other biogeochemical models. In particular, the step size control utilized time step 2 ∆t over the entire simulation for initialization with larger time steps for the NPZ-DOP model.
The tolerance τ 0 ∈ R >0 influenced only the time step's increase for the N, N-DOP and MITgcm-PO4-DOP models. Using smaller tolerances, the possible increase to the maximal time step required more model years for the N, N-DOP and MITgcm-PO4-DOP models. The use of the smallest tolerance always resulted in time step 1 ∆t for the entire spin-up. For the N-DOP and MITgcm-PO4-DOP models, there existed a tolerance for   Table 2) and the N model using the decreasing time steps algorithm (Algorithm 3). Shown are the convergence of the spin-up (left), i.e., the norms of difference (18) between consecutive iterations, and the relative error (29) (right) using different tolerances ε ∈ {0.001, 0.0001}. Furthermore, the figure of the norm of difference (left) contains the convergence towards a steady annual cycle for the reference solution. In this figure, it is mostly covered by that of the decreasing time steps algorithm with tolerance ε = 0.001. which the step size control increased the time step only up to 8 ∆t (Table 7). For the NP-DOP, NPZ-DOP and NPZD-DOP models, the algorithm did not increase the step when tolerance τ 0 = 1.0 was used. Consequently, the choice of a smaller tolerance had no effect (Table 7).

Algorithm 3: Decreasing Time Steps
The decreasing time steps algorithm computed a reasonable approximation of the steady annual cycle with the automatic stepwise reduction of the time step and, thus, reduced the computational effort. For the N model with the parameter vector listed in Table 2, Fig. 5 demonstrates the similar convergence behavior towards a steady annual cycle using the decreasing time steps algorithm as well as the spin-up calculation of the reference solution. Using the decreasing time steps algorithm, the six peaks in the norm of differences (18) pertained to the six reductions of the time step. As a result of the decreased time step, a large concentration change took place in one model year similar to the large changes at the beginning of each spin-up. In particular, the reduction of the relative error resulting from the decrease of the time step is evident in Fig. 5 (right) using the tolerance 0.0001. Although the algorithm applied time step 1 ∆t at the end of the spin-up computation, the spin-up did not perfectly converge against the reference solution. Nevertheless, the algorithm shortened the runtime of the simulation, and finished with an approximation of the steady annual cycle that was much better than the one of the spin-up with constant time step 64 ∆t, i.e., y 10000,64 (cf. Fig. 3). The used tolerance affected only slightly the accuracy of the approximation, but the use of a smaller tolerance resulted in a cost saving because the decreasing time steps algorithm used each time step longer before the time step was reduced. Namely, the application of the decreasing time steps algorithm with tolerance 0.001 resulted in a cost saving of 23% compared to the spin-up computation of the reference solution and in a cost saving of 38% using the tolerance 0.0001. The results for the other biogeochemical models using the reference parameter vectors are parallel to the results shown for the N model.
The decreasing time steps algorithm avoided a divergent spin-up calculation but did not decrease the time step to the minimum for each simulation. We analyzed the decreasing time steps algorithm using the parameter vectors of the Latin hypcercube sample for the different biogeochemical models. For the N, N-DOP and MITgcm-PO4-DOP models, the decreasing time steps algorithm computed an appropriate approximation of the steady annual cycle for the different configurations of the tolerance ε ∈ {0.001, 0.0001} and the number of model years n s ∈ {50, 100, 500} as detailed in Fig. 6 (cf. Fig. 3). More specifically, the accuracy decreased slightly, on the one hand, with a smaller tolerance and, on the other hand, with a larger number of model years n s because in both cases larger time steps were used over a longer period of model years. Figure 7 reflected this behavior in the cost saving. The computational effort increased as a result of a more frequent check of the time step reduction or with a smaller tolerance range because the time step tended to be reduced earlier. For the NP-DOP, NPZ-DOP and NPZD-DOP models, the approximations calculated with the decreasing time steps algorithm were often identical for the different configurations. The algorithm decreased the time step to 1 ∆t only in half of the simulation runs, while this occurred in more than 90% of the simulations using the other three biogeochemical models. In fact, the criterion checking the decrease of the time step was inappropriate if the norm of differences (18) oscillated, cf. [25]. The algorithm, therefore, applied large time steps during the entire spin-up, wherefore the relative error was high (Figs. 6 and 3). In some of these cases, the algorithm temporarily decreased the time step. However, this hardly effected the accuracy of the approximation (Fig. 7). However, the decreasing time steps algorithm automatically reduced the time step if the simulation diverged due to a too large time step. In contrast to the simulations with a low accuracy of the approximation using the NP-DOP, NPZ-DOP and NPZD-DOP models, there were also parameter vectors for which the decreasing time steps algorithm calculated a reasonable approximation with reduced computational costs (Fig. 7).

Conclusions
We presented and evaluated three algorithms that automatically adjust the temporal step size in the spin-up of marine ecosystem models. The setting we used was an offline computation where the ocean transport is realized by transport matrices. We tested the algorithms for six biogeochemical models with different complexity. The first two algorithms are based on a mathematical estimate of the local error taken over a given number of model years. The algorithms are able to both increase and decrease the step size during the spin-up, depending on this error estimate. The first algorithm ignores any negative tracer concentrations, whereas the second adjusts the step size such that negative values are avoided. The third algorithm starts with a coarse step size and reduces it along the convergence of the spin-up. It does not use any mathematical error estimate.
First of all, all algorithms properly worked in the intended way for the investigated model configurations. The spin-up calculations yielded approximations of the steady annual cycle which were in agreement with the respective reference solution obtained with the standard time step. This holds for the considered hierarchy of six biogeochemical models and, for each model, for a reference parameter vector as well as for a Latin hypercube sample of 100 parameter vectors.
We now discuss the behavior of all three algorithm and draw conclusions for them separately. Algorithm 1 automatically adjusts the time step based on the error estimate, but does not take into account any negative tracer values that occur. The error estimation itself needs an additional time-step and, thus, creates an overhead. For the three simpler biogeochemical models, this algorithm chose the biggest feasible step size directly at the beginning of the spin-up. As a consequence, the reduction of the computational effort was very high (up to 95%), since at the upper limit of the step size no further error estimation was performed, and the mentioned overhead did not occur. In contrast, for the three biogeochemical models with higher complexity the algorithm always used the smallest step size. Hence, no performance gain due to bigger time-steps was realized, and, even worse, the error estimation caused maximal overhead. There are some options to vary the behavior of the algorithm. One is the length of the error estimation period. In our computations, this made no qualitative  difference with respect to the chosen step size, but only to the overall performance gain. Two other options that can be easily realized, but were not investigated here, are to estimate the error either only in the beginning or in other selected temporal sub-intervals of the spin-up, but not all the time. Both options will definitely reduce the overhead. However, an adaptive step size control is designed to react on different behavior in different time intervals. Thus, a restriction of the error estimation on certain time pre-defined sub-intervals in the spin-up is somehow arbitrary. It would require some prior knowledge or assumption of the convergence behavior during the spin-up. Another option would be to apply an adaptive time-stepping within one year, and use this choice throughout the whole spin-up. This will also reduce the overhead, since the period for error estimation would be drastically shortened.
Algorithm 2 which additionally strictly avoids negative tracer values did not lead to any performance gain. In all our experiments, it always selected the smallest feasible time-step. Again the error estimation led to a huge overhead. This result shows that these kind of spin-up calculations rely on the acceptance of (small) negative values, even though, in theory, the solutions should stay non-negative for all times. It can be seen that the spin-up results are reasonable even with small negative values occurring at some spatial points.
Algorithm 3 which decreased the time-step during the spin-up depending on its convergence to the steady annually periodic state does not suffer from any computational overhead. It reduces the runtime for the spin-up significantly. The actual reduction depends on the chosen model and its parameters.
As a summary, we can only recommend to use Algorithm 3 directly for an actual spin-up computation. This algorithm is also the easiest to realize. The failing of Algorithm 2 with respect to runtime reduction rises the question how important local negative tracer values are at all for the spin-up runs. This is a question that is beyond the scope of this paper. Finally, from the results for the step size control with Algorithm 1 we deduce that the simpler models do not show different temporal behavior at all during the spin-up, and the more complex do not allow to decrease the temporal resolution. These are results about the models and their behavior in a spin-up: For the simpler models, this means that their spin-up convergence is rather uniform. However, for other existing far more complex models (see, e.g., [10]) the question remains if these models also show this uniform spin-up convergence. For such kind of models, Algorithm 1 might be a useful tool to study the convergence behavior during the spin-up.