Earth System Model Parameter Adjustment Using a Green’s Functions Approach

We demonstrate the practicality and effectiveness of using a Green’s functions estimation approach for adjusting uncertain parameters in an Earth System Model (ESM). This estimation approach had previously been applied to an intermediate-complexity climate model and to individual ESM components, e.g., ocean, sea-ice, or carbon-cycle components. Here, the Green’s functions approach is applied to a state-of-the-art ESM that comprises a global atmosphere-land configuration of the Goddard Earth Observing System (GEOS) coupled to an ocean and sea-ice configuration of the Massachusetts 5 Institute of Technology general circulation model (MITgcm). Horizontal grid spacing is approximately 110 km for GEOS and 37–110 km for MITgcm. In addition to the reference GEOS-MITgcm simulation, we carried out a series of model sensitivity experiments, in which 20 uncertain parameters are perturbed. These “control” parameters can be used to adjust sea-ice, microphysics, turbulence, radiation, and surface schemes in the coupled simulation. We defined eight observational targets: sea-ice fraction, net surface shortwave radiation, downward longwave radiation, near-surface temperature, sea surface temperature, sea 10 surface salinity, and ocean temperature and salinity at 300 m. We applied the Green’s functions approach to optimize the values of the 20 control parameters so as to minimize a weighted least-squares distance between the model and the eight observational targets. The new experiment with the optimized parameters resulted in a total cost reduction of 9% relative to a simulation that had already been adjusted using other methods. The optimized experiment attained a balanced cost reduction over most of the observational targets. We also report on results from a set of sensitivity experiments that are not used in the final opti15 mized simulation but helped explore options and guided the optimization process. These experiments include an assessment of sensitivity to the number of control parameters and to the selection of observational targets and weights in the cost function. Based on these sensitivity experiments, we selected a specific definition for the cost function. The sensitivity experiments also revealed a decreasing overall cost as the number of control variables was increased. In summary, we recommend using the Green’s functions estimation approach as an additional fine-tuning step in the model development process. The method is 20 not a replacement for modelers’ experience in choosing and adjusting sensitive model parameters. Instead, it is an additional practical and effective tool for carrying out final adjustments of uncertain ESM parameters.


Introduction
Earth System Models (ESMs) include various parameters that govern the representation of unresolved, unrepresented, or under-observed processes in the models. The most sensitive parameters are typically adjusted during the last step of model 25 development relative to observational targets. Currently, there is no agreed-upon methodology to adjust model parameters. As an illustration of the range of approaches and observational targets, Schmidt et al. (2017) describe the vast range of tuning strategies used for ESM development in U.S. modeling centers. Parameter adjustment aims to improve the representation of various processes or key state fields in the model relative to a predetermined set of observational targets (Mauritsen et al., 2012;Hourdin et al., 2017). 30 ESM tuning is often done in a heuristic, trial-and-error way, wherein one or more parameters are perturbed to new values based on diagnosing the causes of systematic error with the aim of improving the model's overall behavior (e.g., Watanabe et al., 2010;Yukimoto et al., 2012;Mauritsen et al., 2019;Voldoire et al., 2019;Golaz et al., 2019;Sellar et al., 2019;Danabasoglu et al., 2020). The observational targets can include, for example, radiation balance at the top of the atmosphere, surface temperature, sea ice variability, tropospheric wind, or even a derived target such as a "cloud object" (Posselt et al., 2015). In this 35 approach, a typical tuning exercise comprises several iterations of suites of sensitivity experiments with one or more perturbed parameters. After each iteration, a comparison is made against a set of observations, and another set of sensitivity experiments targeting either the same or another set of parameters and observations is conducted based on the results. This tuning exercise is highly dependent on the experience and expertise of model developers and on the application for which the model will be used. 40 There are three main drawbacks for heuristic optimization approaches. The first is that this method necessitates an analysis step between each sensitivity sweep, a time-consuming process. The second drawback is the large number of needed simulations, since a new set of experiments is required for iteration of the process. The third drawback is the interdependence of each parameter's effect on the resulting ESM simulation. For example, both sea ice albedo and precipitation efficiency could affect the globally-averaged radiation balance at the top of the atmosphere. The adjustment of one parameter or set of parameters at a 45 time will be suboptimal because estimates of empirical parameters depend on each other and on model configuration, boundary conditions, etc.
Given unrestricted computer resources, one could randomly explore model parameter space exhaustively and choose the simulation that best fits all available observations. In practice, this type of exhaustive parameter exploration is not feasible and we need tractable methods that combine objective methodologies with the model developer's experience to arrive at an 50 optimized choice of model parameters.
To mitigate the drawbacks and computational cost of above approaches, a number of more objective parameter calibration methods have been proposed for systematically choosing or "fine-tuning" the final set of optimized parameters. These methods include the use of Latin Hypercube sampling techniques (Posselt et al., 2015), the downhill simplex method (Zhang et al., 2015), the multiple very fast simulated annealing (Zou et al., 2014), and Gauss-Newton line search algorithms (Tett et al., 55 2013(Tett et al., 55 , 2017Roach et al., 2018). The above estimation methods require a large suite of experiments, which are computationally expensive. Therefore, the sensitivity experiments are typically carried out with very low resolution, atmosphere-only, or even single-column model configurations. The direct applicability of these methods for tuning an ESM may consequently be limited.
Another alternative to heuristic estimation approaches is the adjoint method. Using the adjoint method, one can adjust many parameters simultaneously. A successful example is provided by the Estimating the Circulation and Climate of the Ocean 60 (ECCO) consortium, which used the adjoint method to adjust Massachusetts Institute of Technology general circulation model (MITgcm) ocean model parameters, initial conditions, and boundary conditions (Forget et al., 2015). Although the adjoint method provides robust results, developing a model's adjoint is time-consuming, even when using an automatic differentiation software tool. Additionally, the adjoint method is computationally expensive because many forward-adjoint ESM simulations would be required to optimize a non-linear model.
The model used for this study is the GEOS-MITgcm coupled model. We briefly describe here the particular configurations of these two models as they are used in our study. The GEOS AGCM's dynamical core and suite of physical parameterizations, along with the land model and aerosol model are described in Molod et al. (2020) and in the references therein. The turbulent surface layer parameterization, relevant for this study, is a modified version of the parameterization documented in Helfand 95 and Schubert (1995). The wind stress and surface roughness model was modified by the updates of Garfinkel et al. (2011) for a mid-range of wind speeds and further modified by the updates of Molod et al. (2013) for high winds. The GEOS AGCM's cubed-sphere grid was configured to run with a nominal horizontal grid spacing of 100 km. The vertical grid is a hybrid sigma-pressure grid with 72 levels. The time step of integration was set to 450 seconds.
MITgcm has a finite-volume dynamical core (Marshall et al., 1997). It has a nonlinear free-surface and real freshwater flux 100 (Adcroft and Campin, 2004) and a nonlocal K-profile parameterization scheme for mixing (Large et al., 1994). The MITgcm horizontal grid is the so-called "Lat-Lon-Cap" (Forget et al., 2015), and was configured to run with a nominal grid spacing of 100 km. The vertical grid is the z * height coordinates (Adcroft and Campin, 2004), and it has 50 vertical levels. The time step of integration was set to 450 seconds, similar to the atmospheric model.
The GEOS-MITgcm atmosphere-ocean interface includes a skin layer (Price et al., 1978), configured for the simulations 105 described here to impose a one-day time scale on the interaction, and the communication between the ocean and atmosphere is updated at every atmospheric model time step (450 seconds). The sea-ice dynamics model is provided by MITgcm and the sea ice thermodynamics model is from CICE 4.1 (Hunke and Lipscomb, 2010).

Mathematical Formalism and Methodology
Algebraically, the ESM described in previous section can be expressed as a set of rules for time stepping a state vector, where x(t i ) is a discretized representation of the Earth System at time t i and function M i represents ESM time-stepping rules for advancing the state vector from time t i to t i+1 . The discretized evolution of the true Earth System, x t (t i ), is assumed to differ from that of the numerical model by a vector of stochastic perturbations η so that where ε represents a vector noise process assumed to have zero mean and covariance matrix R. In addition to measurement errors, vector ε includes all model and estimation "representation" errors, i.e., model-data residuals that cannot be represented 120 by M i and η.
To fit the ESM to the available observations, we aim to minimize an objective cost function which penalizes a weighted least-squares distance between optimized and prior model parameters (η T Q −1 η) and a weighted least-squares distance between simulation and observations (ε T R −1 ε), where T is the transpose operator and −1 indicates matrix inversion. It can be shown that the minimization of Equation (4) will produce the maximum likelihood solution for Gaussian inverse problems Kalnay (2002). For the Green's Functions approach, we adjust a small, carefully-chosen set of parameters η. To do this, we combine Eqs.
(2) and (3) to obtain: where G represents a function the combination of observation operator H with ESM time-stepping rules M i . We assume that 130 (5) can be linearized about a reference ESM simulation so that: where 0 is a null vector, G(0) represents the reference ESM simulation sampled by measurement function H, vector y d is the model-data difference between observations and the reference simulation, and kernel matrix G can be thought of as a first order multivariate Taylor expansion of G around η = 0.
The assumption here is that for small enough perturbation, the model 135 response is linear. Kernel matrix G can be derived by conducting an ESM perturbation experiment for each element of vector η. Specifically, the j-th column of matrix G is where e j is a perturbation vector filled with zeros except for element j, which is perturbed by scalar value e j . The minimization of (4) given (5) is a discrete linear inverse problem with solution: where is the uncertainty covariance matrix for the optimized estimate η a . Mathematically, the Green's function approach is equivalent to any linear, Gaussian inverse problem, for example, the adjoint method.

145
In practice, the Green's functions estimation method can be divided into six steps: 1. Run a reference ESM simulation using default values for model-error parameters, that is, η = 0. The output from this experiment, sampled by observation operator H, provides G(0) in (5)  linear combination and a projected cost under the assumption of linearity.
6. Run a new simulation using optimized parameters η a . This optimized simulation can be compared with the optimized linear combination of sensitivity experiments in order to evaluate the linearization assumption.
Note that steps the choice of the observational targets and the definition of the covariance matrices (steps 2 and 3) can be later revised without the need to rerun a new set of experiments (step 4).

Proof of Concept Optimization
In this section, we illustrate the application of the six steps listed above.

Reference Experiment
The 10-year long reference experiment was configured with reference values of a set of parameters that we aimed to optimize.
The model was configured to run in "perpetual year" mode, meaning that the external boundary conditions (solar insolation, 165 greenhouse gas amount and aerosol emissions ) are kept to those of 2000. The "perpetual year" mode was chosen to optimize the model's equilibrium state relative to our climatology observational targets.

Observational targets
The motivation of this study was to set up a simple, systematic, reproducible, extensible, and efficient framework for improving the climatology of the new coupled model now and as it undergoes development in the future. To start, we chose eight key 170 observational targets (Table 1)

180
In this study, Q −1 is set to zero, that is, the cost function does not explicitly include a priori knowledge about the range of possible parameter values. Some of the parameters, however, have intrinsic physical bounds for the range of values. For these parameters, the resulting optimized parameters were checked and all were found to lie within the range of physically-acceptable values. The error covariance matrix R is chosen to be diagonal. Because the estimation problem is highly overdeterminedthe number of parameters to adjust is O(10) while the number of observations is O(10 6 ) -we expect that the assumption 185 that R is diagonal will have minimal impact on the optimized estimate η a . In general, the addition of non-diagonal elements of R would indicate linear dependence among observational errors, reducing the effective degrees of freedom. Since, in our case, the number of observations is far larger than the number of optimized parameters, we would still retain a very large number of independent observations as compared to optimized parameters. The variance was calculated separately for each observational target and defined to be the global, area-weighted, mean-squared difference between the climatology of the 190 reference experiment and observations, that is: where the indices v, i, j, and s represent observation variable, longitude, latitude, and season (winter, spring, summer, and fall), the hat signˆrepresents climatological seasonal-mean values, and w i,j is an area weight with i,j w i,j = 1/4, giving each season a 1/4 of the weight. The prior error variance of each variable is, therefore, a constant value. A consequence of 195 our choice of the covariance matrix and the cost function is that the cost of the reference experiment is equal to one for each variable. For eight observational targets, the total cost of the reference experiment is eight. This definition gives equal weight to each of the eight observational targets and was found to provide a balanced solution in terms of each observational target's overall influence. The cost of the optimized experiment is expected to be smaller than one for each variable. In addition to the reference experiment, 20 10-year long sensitivity experiments were performed using the exact model configuration as the reference experiment but perturbing one of the uncertain parameters each time. The length of the simulations needed for the Green's Function optimization is related to the application of the model being optimized. For climate projec-tion applications, where the long-term equilibrium of the model is paramount, it seems clear that longer simulations would be warranted in order to allow for model spin up time. The primary application of the coupled model presented here is seasonal 205 to decadal prediction, and so the optimization is done to include (ideally to minimize) the initial model drift.

Sensitivity Experiments
The reference and perturbed parameter values, and a short description of the parameters is listed in Table 2. We used the 20 sensitivity experiments to estimate the model's sensitivity to each of the perturbed parameters based on the Green's functions formalism. The perturbed parameter values where chosen based on expert advise. For some of the parameters the value was chosen based on physically realistic values from past experience. In some of the cases, the perturbed value was chosen 210 arbitrarily due to a large parameter uncertainty. Note that there is no risk of overfitting due to the large number of observational targets (eight observational targets, two-dimensional fields, four seasons -≈300,000 targets). Methods to penalize a large number of parameters are therefore not needed in this case.

Optimized Parameters
The optimized parameters are given in the last column of    Table 3) shows that the total cost reduction of the optimized experiment (shown on the right) relative to the reference experiment (shown on the left) is 9% (similar to the cost reduction found in Zhang et al. (2015)). Figure 1 (and Table 3) also shows the projected cost, which is the cost one would get if the model response to parameter perturbations were entirely linear. Most of the cost reduction seen in the projected (second from the right) experiment is realized in the optimized 230 experiment. The total cost of each of the 20 sensitivity experiments is larger than the cost of the simulation with the optimized parameters. Figure 2 shows the normalized cost for each of the eight observational targets in the same order as in Figure 1.  longwave radiation. The longwave radiation had a cost increase of 2%. The more considerable increased cost of 16% percent in the surface temperature is related to our choice of the observational targets. Five out of the eight observational targets are ocean-only variables, only one is land-only, and two cover both land and ocean. This asymmetry between land and ocean constraints seems to have pulled the parameters to values that benefit the ocean more than the land. This proof of concept 240 illustrated in this study, although clearly demonstrates the overall cost reduction with this method, also illustrates the need to choose the desired observational targets carefully. The cost reduction was also found to be rather homogeneous across many regions around the globe. Figure 3 shows that, for net shortwave radiation, the total cost reduction was 20% relative to the reference experiment. Most of the oceans exhibit a large reduction in the cost, particularly over the Antarctic Circumpolar Current (ACC) and in the stratocumulus regions in the The actual cost reduction of the various variables seen in the optimized experiment is also generally consistent with the projected cost. Figure 4 shows the projected cost reduction for the net shortwave radiation, which is in broad agreement with the actual cost reduction of the optimized experiment seen in figure 3c. In this case, the cost of the optimized experiment was 250 even smaller than the projected cost. Overall agreement between the projected cost and optimized experiment cost suggests that the linearity assumption that underlies the Green's functions methodology holds sufficiently well for the approach to be useful.
The surface temperature, which exhibits an overall 16% increase in cost ( Figure 5), also shows some large regions of cost reduction over North America and Greenland. However, the cost in north Asia and Europe is increased. In terms of the 255 cost's seasonality, it was found that DJF and MAM had the highest cost increase while JJA months had the smallest overall cost reductions. These results suggest that the misrepresentation of land-ice processes in the new and the old configurations is responsible for the high cost. Surface temperature is the variable that had the least projected cost reduction. So, it was expected to exhibit the worst performance in terms of the optimized experiment cost. Including more land-related observational targets such as snow fraction and soil moisture could improve the cost of surface temperature. This inclusion can assist in cost increase, but here the increase is 2%. In terms of the other variables, the cost of sea-ice fraction was reduced by 8.2%, and the four ocean-only variables showed cost reduction between 10% to 25%, indicating an improved representation of the ocean circulation.

265
The parameter estimation experiment presented in the previous section describes our current best practice, which was guided by a series of sensitivity experiments with different configurations where different variations of the methodology were tested.
Below we discuss sensitivity to the choice of (1) prior covariance matrices, (2) observational targets, and (3)

Sensitivity to choice of prior covariance matrices
In general, the error covariance matrices Q and R are not known and difficult to estimate. For the Green's function approach, however, the number of observations is generally much larger than the number of parameters being estimated, which simplifies the selection of Q and R. First, the small number of control parameters limits the solution's degrees of freedom; therefore 275 the choice of Q and R, if they are reasonable, is not expected to change the solution η a much -but it will impact posterior uncertainty P. Second, the data kernel matrix G is small enough to be defined explicitly; therefore many interesting properties of the solution can be derived and evaluated. Third, the solution of (8) and (9), once the kernel matrix G has been derived, is trivial; therefore it is possible to test the impact of particular choices of Q and R, as is done next.
The sensitivity to the choice of covariance matrices Q and R in Equation (4) is evaluated using a pared-down configuration, 280 perturbing only three parameters. The pared-down configuration is used to simplify the discussion and better illustrate the response of the optimized solution to different covariance matrices. The perturbed parameters were AUTscale, LTS_LOW, and two sea ice parameters (ALBICEV and ALBICEI) perturbed together. The output from the methodology for this choice is a set of scaling factors for the parameters. When more than one parameter is perturbed in one experiment, the same scaling factor is applied to all the experiment parameters. In this configuration, unlike in our best practice configuration, the observational 285 targets were multi-year means, and the cost was not computed separately for each season. The observational targets are also similar to the best practice configuration except that the 500mb height was used in place of the surface temperature over land, and the net radiation replaced the two radiation components. We tested four different options for the prior error variance, the diagonal elements of matrix R: where, v, i, j, s, t are the variable, longitude, latitude, season, time indexes. The tilde symbol (˜) represents seasonal means and N is the number of years. The first option (a) uses the area-weighted variance of the observations instead of the mean-squared 295 misfits used in Eq. 10. The second option (b) is the spatially and seasonally variable, time variance of the observations, the third option (c) is the square of the maximal absolute difference between observations and the reference experiment, and the fourth option (d) is an area weighted variance of the difference between the observations and the reference experiment. All four estimates of R are calculated separately for each variable. The four estimates of R aim to represent a range of plausible estimates for the actual error variance and allow us to explore the effect on the resulting parameters and cost. 300 We calculated optimized parameters η a for each of the four definitions of R assuming that Q = ∞, that is, no prior information about the control parameters, and again assuming that Q = I, the identity matrix, meaning that we have the same amount of confidence about the prior values of each control parameter. These two estimates of Q represent two extreme situations. Table 4 shows the eight different combinations of the parameters based on the four options for the observations covariance and two different parameters governing the covariance matrices. In general, all configurations generate realistic values in terms 305 of the order of magnitude, but parameters differ in detail. The Q = I case tends to keep parameters closer to their initial value (relative to Q = ∞) as expected. Table 5 shows the cost associated with each of the eight experiments. Looking at the Q = ∞ cases, defining the cost based on Eq. 11a reduce the cost by 7.7% but the relative importance of each variable varied substantially. For example, 500mb height cost was very small for the reference (0.024) and the fact that the optimized experiment increase the cost by 33.7% did 310 not affect the total cost. Also, the sea ice contribution to the cost was much larger than the temperature and salinity at 300m below the sea surface. Option 11b was found to provide only marginal cost reduction. Option 11c demonstrated the largest cost reduction (the smallest ratio in the J_total column). However, most of the projected cost reduction was restricted to the 500mb height target, while the cost of most of the other variables was not reduced, particularly in the case of the ocean targets. This is in contrast with option 11d, for which the cost reduction was spread more uniformly across the observational targets. All of 315 these considerations led to our choice of option 11d as the best practice. The Q = I was not chosen because the range of the parameters is in itself often quite uncertain and because optimized parameter values were found to remain reasonable without the need for adding this additional constraint to the cost function. This choice could certainly be worth revisiting in a future study.
After running a new model simulation with the optimized parameters based on Eq. 11d and Q = ∞ case, the total cost 320 reduction was found to be 6%, 2.9% less than the projected value, but still not small (Table 5). The 500mb height contributes to most of the reduction, and the seaice cost was increased by about 25%. A possible explanation for the 500mb target being the main contributor to the cost reduction is that there were a small number of grid points in which the error was reduced substantially, particularly at high latitudes (not shown). This result suggests an experiment in which a spatially-dependant weight is used to reduce the magnitude of the outliers. The current cost implementation has the dependency on the model 325 squared-error. This dependency favors the correction of a small number of points with large errors relative to a large number of points with small errors. Choosing absolute error instead of squared-error would have changed the tendency to favor the correction of a small number grid points with large errors. Nevertheless, the simpler experiment configuration was a successful choice to explore the sensitivity to the covariance matrix designs.  Here, we evaluate the sensitivity of the cost to the number of optimized parameters. We calculated the projected cost reduction for four cases which used a subset of the first 5, 10, 15 and all 20 parameters from Table 2. We focus here on an optimization based on the multi-year means. An increased number of parameters is expected to translate into a reduced cost as we fit the model to more parameters. In the following results we decided to remove the 500mb height from the cost function in order to spread the cost across more variables. Generally, there is nothing wrong with having most of the cost reduced in one of 335 observational targets, but we decided here to seek a more balanced cost reduction. An additional change here and in the following sections is the separation of the net surface radiation target into net shortwave and downward longwave radiation, as they are independent and both have SRB observational counterparts. The choice of downward longwave rather than the net was made due to the direct dependence of upward longwave radiation on the SST (which is already a constraint).
We tested the projected cost reduction as a function of the optimized parameters' number (figure 6). We found that the 340 projected cost was reduced from 6% to 25% when going from 5 optimized parameters to 20. The cost reduction was gradual showing and additional reduction with the increase of the optimized parameters. The largest cost reduction was found between 5 to 10 optimized parameters.

Annual-Mean versus Seasonal-Mean Data Constraints
Here, we used nine sensitivity experiments, with nine optimized parameters, to look at the sensitivity of the optimized pa-345 rameters to seasonality. Two sets of parameters were calculated -based on multi-year means (annual, Table 6) and based on multi-year seasonal means (seasonal, Table 6). The calculation of the optimized parameters comes after the simulations were performed, so the input simulations are the same for both options. The only difference between the annual and seasonal suite is the configuration of the cost function provided to the Greens' functions methodology.
The most noteworthy difference between the annual and seasonal results was found in the Charnock parameter that controls 350 surface winds. Although the difference between parameters is smaller than 10%, the projected cost differences are large, with the annually-based optimization reducing the cost by about 15% and the seasonally-based optimization reducing the cost by about 9%. This result is expected, as seasonally-based data have more variability and we use the same number of parameters to optimize model results to larger number of observations.
Optimizing for seasonal observational targets is expected to be essential for improving seasonal variability. Therefore, de-355 spite the less efficient cost reduction for the seasonal targets we decided to continue to focus on experiments with seasonal observational targets. For these two sets of experiments, we did not choose to rerun the model with optimized parameters, but instead decided to increase the number of optimized parameters by performing additional simulations.

Length of Optimization Period
This sensitivity test was done with the 20-parameter set of simulations. The appropriate length of the sensitivity experiments 360 was chosen based on the intended application of our tuned model for seasonal and decadal prediction. It seems clear that for a different application of the model, say climate projection, the length of the sensitivity experiments, optimization and  optimizing the latter period only is that 2005-2010 had fewer than two Elnino cycles, and it was too short for the calculation of the model climatology. The 10-year experiment had three full Elnino cycles, and its climatology represents the model climatology more adequately. The 10-year optimization, therefore, is found to enable better parameter optimization.

Summary and Concluding Remarks
This study demonstrates the applicability of the Green's functions approach, introduced in Menemenlis et al. (2005) for a 375 forced ocean simulation, to the adjustment of uncertain parameters for a coupled ocean-atmosphere simulation. In this proof-ofconcept study, we computed the response of the coupled ocean-atmosphere simulation to the perturbation of the 20 parameters listed in Table 1. These 20 perturbation experiments were subsequently used to optimize the values of the 20 parameters.
The observational targets were long term seasonal-means of the eight fields listed in Table 2. The optimization increased the explained variance by up to 25% (for SST in Table 3).
The Green's functions approach assumes that the response of the model to small perturbations is approximately linear. This assumption does not rigorously hold for atmospheric, oceanic or coupled-ocean-atmosphere models because of non-linear weather and climate phenomena such as synoptic storms and El Niño-Southern Oscillation (ENSO) variability. Therefore, a key consideration is that the perturbation experiments be sufficiently long to average out non-linear weather and lower frequency phenomena. This study shows that 10-year long simulations can be sufficient to enable successful application of the 385 Green's functions approach.
The cost reduction was spread across most of the observational targets, though our configuration of observational targets may have been responsible for the cost increase in two amongst them -the targets associated with land areas that are underrepresented in our cost function. In the end, The choice of the observational targets should reflect the objective of the study and, probably, there is no single set of observational targets that is adequate for all applications a priori. Our proof-of-concept study 390 and the experiments that show the impact of different choices in the details readily reflect the trade-offs among data constraints that can be important considerations.
Increasing the number of optimized parameters reduced the overall cost which seems to be quasi-linear with the number of parameters, but further investigation into this is required. Increasing the temporal resolution of the observational targets from annual to seasonal increased the cost further -a result of the increasing number of the observational targets. Adding 395 seasonality in parameter adjustments is viewed as one of the next logical steps in this regard. We also tested four different methods to calculate the error covariance matrix, and found that using the error of the reference experiment produced the most balanced results in terms of accounting for all observational targets.
This study shows that the Green's functions methodology can benefit the earth system modeling community by providing a more structured yet practical method to tune the complex global models. The Green's functions methodology has several Identifying the most important uncertain parameters for applying Green's functions methodology still requires close familiarity with models and there is no replacement for the experience of the modelers when using this methodology. Therefore it is viewed as a valuable framework to further leverage modelers' expertise to fine-tune model parameters, standardize practices 410 across the various modelling centers, and improve model products beyond the present state of art.
Code and data availability. Data leading to this publication was made available as part of a Zenodo repository at: https://doi.org/10.5281/zenodo.5507631