Computer model calibration based on image warping metrics: an application for sea ice deformation

Arctic sea ice plays an important role in the global climate. Sea ice models governed by physical equations have been used to simulate the state of the ice including characteristics such as ice thickness, concentration, and motion. More recent models also attempt to capture features such as fractures or leads in the ice. These simulated features can be partially misaligned or misshapen when compared to observational data, whether due to numerical approximation or incomplete physics. In order to make realistic forecasts and improve understanding of the underlying processes, it is necessary to calibrate the numerical model to field data. Traditional calibration methods based on generalized least-square metrics are flawed for linear features such as sea ice cracks. We develop a statistical emulation and calibration framework that accounts for feature misalignment and misshapenness, which involves optimally aligning model output with observed features using cutting edge image registration techniques. This work can also have application to other physical models which produce coherent structures.

Sea ice is frozen sea water that insulates the warmer ocean from the colder atmosphere in winter and provides a partial barrier to heat, moisture, and momentum transfer between the ocean and atmosphere, making it an important component in the Earth's energy balance. The ice cover waxes and wanes seasonally and at its maximum extent covers around 7% of the Earth's surface and close to 12% of the Global Ocean Weeks [2010]. Cracks in the ice cover, called leads, occupy 1-2% of the ice cover in winter but account for half of the heat flux from the ocean to the atmosphere Badgley [1961].
A constitutive model, the elastic-decohesive model [Schreyer et al., 2006] has been developed to describe crack formation in sea ice. The important aspect of this model is that the existence as well as the orientation and the opening or closing of a crack is directly represented.
Parameters defining the material properties for the elastic-decohesive model were often chosen based on experiments [Schulson, 2004] and some in situ observations [Coon et al., 1998]. There has never been a more detailed study of the appropriateness of these parameters and whether other choices would better fit observations.
With new high resolution satellite data now available, unprecedented opportunities to improve our understanding of the interaction between dynamics and thermodynamics of sea ice, and to appraise new models are now possible. With this imagery, we are able to resolve long linear features in the pack ice that evolve with time that can generally be associated with fractures or leads. These leads are derived from the small-scale, time-varying deformation of the ice cover observed using the RADARSAT Geophysical Processing System (RGPS) [Kwok, 2001].
Thus, a fresh assessment, through calibration of the model parameters, is appropriate. The focus of this paper is to develop metrics suitable for comparing simulation and observations of these linear features in sea ice, and to use these metrics to calibrate the elastic decohesive model. Calibration is the traditional statistical method to infer unknown input parameters for a model of a physical process from the corresponding observed output. The Bayesian framework O'Hagan, 2001, Higdon et al., 2004] is a natural choice for calibration where the inference is based on the posterior distribution of the model parameters given the observed outputs. Here, the likelihood term in the posterior will include the non-linear physical forward model, via, the elastic-decohesive model, which is computationally expensive to solve numerically. Moreover, the likelihood needs to be evaluated at each iteration of the Markov Chain Monte Carlo (MCMC) sampling scheme, thus making the Bayesian calibration method computationally very challenging. Alternatively, a statistical emulator can be constructed using the simulation runs of the elastic-decohesive model at some suitably chosen design points.
The emulator can then act as a surrogate, providing fast predictions of the simulation output at unsampled input values, with corresponding measures of uncertainty O'Hagan, 2002, Conti et al., 2009]. In practice, a Bayesian modularization approach [Bayarri et al., 2007, Liu et al., 2009] is often used, in which a Gaussian process (GP) based emulator is fit only once, then the fitted emulator replaces the physical model inside the likelihood for Bayesian calibration. We adapt this in our proposed method.
There are two major challenges in applying the traditional emulation-calibration framework to linear features such as sea ice cracks. (1) In the traditional framework, the output from the elastic-decohesive model would be modeled as a realization from a GP; however, emulating the response surface with linear features using a GP is inadequate. Although the GP model is attractive for emulation because of its flexibility to fit a large class of response surfaces, it is not suitable for emulating a surface with sharp gradients and discontinuities. For model outputs with linear features, a simple GP emulator would tend to smooth out those features. (2) Generalized least-square metrics are implicitly used in the likelihood functions when fitting a GP emulator and carrying out calibration; however, such metrics are also flawed when used in the presence of sharp gradients and discontinuities. Because the linear features occupy only a small, usually lower dimensional, region of the domain, models that completely fail to predict these important features may actually score better by these metrics than those that simulate the dynamically important features albeit at a somewhat different location or orientation than observations indicate [Mass et al., 2002]. We propose a novel emulation-calibration framework, which utilizes cutting-edge image registration technique, to circumvent the above challenges in dealing with linear sea ice features. It involves optimally aligning high-dimensional model outputs to the observed features, then performing emulation and calibration based on the metrics provided from image registration.
Image registration, finding a correspondence between pixels across images, is a standard problem in image processing (see [Goshtasby, 2012] for a review). Recently, the concept of warping images for registration using metrics has become increasingly popular given the properties they provide (Beg et al. [2005], Tagare et al. [2009], Vercauteren et al. [2009], Xie and Srivastava [2016]). Image registration is usually posed as a variational problem where the objective function between two images f 1 and f 2 is of the form where | · | is the Euclidean norm, γ is the deformation, R is a regularization to ensure smooth deformations, and Γ is the space of relevant deformations. Recently, this metric has been applied to model emulation and calibration [Kleiber et al., 2014]. However, the 'metric' used in Kleiber's approach is not a proper metric in the image space under warping and can fail to properly calibrate when observation and simulation images are not properly aligned.
The underlying problem with the above formulation is that it does not provide a proper metric that is inverse consistent. Here we propose a new method using the concepts of phase and amplitude of an image. The power of this method is that we have two metrics, first proposed by Xie and Srivastava [2016], which are proper distances in the space of images. One, the phase energy, gives a measure of how much warping is required to register geometric features. The other, the amplitude energy, gives a measure of the pixel intensity differences, if any, after an optimal warping has been applied. The combination of the measures allows one to determine how well a model is reproducing geometric features or coherent structures. A low amplitude energy, for instance, implies there is a smooth deformation from one image to the other, meaning a geometric feature is likely captured but perhaps misaligned. The magnitude of the phase energy gives a measure of how misaligned the feature is. In contrast, a high amplitude energy implies the geometric features of the two images are inconsistent and are not captured by the model. Using these measures in our calibration scheme allows the emulator to better distinguish how different a model prediction truly is from an observation.
We propose to use these two informative metrics to train our emulator for use in calibration. This allows us to use fewer parameter settings for expensive model runs while providing a unique measure of model performance. The advantage of our approach is two-fold: (1) the alignment avoids penalizing predicted features that are misaligned and misshapen in calibration, and (2) emulation-calibration for scalar metrics is straightforward so that standard software can be applied directly. The rest of this paper is organized as follows, Section 2 introduces the RGPS deformation data and the sea ice model that describes the physics of sea ice deformation. Section 3 introduces image warping and the emulation-calibration framework based on the metrics provided by the warp. We illustrate the usefulness of our approach via a toy example and a perfect model experiment based on the sea ice model in Section 4, and we apply the proposed method to the RGPS data and present the results in Section 5. Some final remarks on the proposed approach are included in Section 6.

Data
The fundamental quantity produced by RGPS is trajectories of marked sea-ice points. At an initial time, a set of points forming a square grid is located in a synthetic aperture radar (SAR) image. In images resulting from subsequent satellite passes, the original points are found again using area-based and feature-based tracking. The time separation between repeat observations is variable and is based on available coverage. The time interval between successive images is called the timestep. This procedure provides the trajectory of each point as these points move with the ice cover. Since the same set of points is tracked over a ice season, RGPS provides a densely sampled Lagrangian picture of the motion, similar to what would be obtained by drifting buoys. Secondary procedures in RGPS derive estimates of ice deformation. If the set of points in the original configuration is viewed as vertices of grid cells then the motion of the points determines the deformation of those cells. With this interpretation, grid quantities can be approximated that help provide a picture of the rate and type of ice deformation.
The RGPS deformation product consists of cell area changes and spatial derivatives computed using the ice displacements at the cell vertices. Each cell has a unique identifier, the time of its creation, the map coordinates of the cell center, the displacement of the cell center between this and the last observation, the timestep and the derivatives of displacement increments. Since leads are formed by cracks in the ice, the displacement associated with a lead is discontinuous. The article by Coon et al. [2007] introduced the idea of separating the displacement into the sum of a continuous and discontinuous part. A kinematic analysis of RGPS deformations was performed assuming that the predominant term was the jump discontinuity in displacement. This procedure fit a jump in displacement, or a crack, in each RGPS cell that could best account for the observed deformation. This fitting procedure provides the orientation and amount of normal and shear opening of the crack. Details on this fitting processes can be found in . Figure 5(b) shows an example of applying the procedure to a region of the Beaufort Sea to plot the magnitude of the predicted jump in displacement.
Observations and simulations of a region of the Beaufort Sea between Banks Island and Point Barrow during the period from Feb. 23 -Mar 11, 2004, are analyzed in this paper. During this period, daily observations are available on a 10 km square grid. Specifically, the domain represents a crack with a normal to the crack given by n = (cos α(s), sin α(s)) and a tangent to the crack given by t = (− sin α(s), cos α(s)) in the plane of the ice. Thus, the orientation of the crack is determined by α, the angle that the normal makes to the s 1 -axis. The size of the crack is obtained from the normal and tangential components of the jump vector, u n and u t , respectively. We threshold the jump in displacement and drop fractures smaller than 0.4 km.
Simulations using the material-point method (MPM) for this problem are described in Sulsky and  and output the same jump vector, [u], as that which can be derived from the RGPS data. To summarize, the simulations are performed on the same 10 km grid as the observations, initially with four material points per cell. Land is modeled with rigid, immobile material points. Boundaries of the rectangular region intersecting ice are given a prescribed velocity in agreement with the RGPS observations. At the start of the calculation the ice is at rest with a uniform average thickness. The thickness distribution producing the average thickness is assumed to be parabolic with a maximum at h = 3 m and equaling zero at h = 0 m and h = 6 m. Six hour wind fields from the National Center for Environmental Protection and National Center for Atmospheric Research (NCEP/NCAR) reanalysis [Kalnay et al., 1996] are used to determine the wind drag. Additionally, six hour air temperature, specific humidity, long wave flux, and short wave flux from the NCEP/NCAR reanalysis are used for the ice column temperature calculation. Details of the thermodynamic model are also in Sulsky and Peterson [2011]. The ocean currents are updated daily using the output from an ocean model (MITgcm, Marshall et al. [1997]) run independently from the ice simulation. Nominal values of sea surface temperature and salinity of −1.8 • C and 32 ppt are used throughout the domain for the calculation.
The elastic-decohesive constitutive model used in the simulations is based on two elastic parameters, Young's modulus, E, and Poisson's ratio (ν = 0.36), and five independent decohesive parameters. There are three strength parameters: tensile failure strength (τ nf ), shear failure strength (τ sf ), and compressive failure strength (f c ); plus the shear magnification factor set to 4, and the opening parameter, u 0 , set to 0.4 km. For the calibration, τ nf is chosen uniformly in the range 10-50 kPa and the shear strength in the range 10-150 kPa. These choices allow the shear strength to be both larger and smaller than the tensile strength.The compressive strength is set to 5τ nf and the Young's modulus to E = 10f c . Keeping these ratios means that u 0 would not need to be changed to maintain stability of the algorithm and avoid snapback. Note that these strength ranges are smaller than what has been seen in experiments (eg. Timco and Weeks [2010]), but are larger than what is typical for large scale simulations. The fact that sea ice on a large scale might be effectively weaker could be attributed to features such thermal bending fractures that appear on a scale of about 1 km, or other impurities and weaknesses that would not show up on a laboratory scale, but are smaller than typical simulation scales.
Fractures are initialized in the simulations using the kinematic analysis of the RGPS data over one day (23 February). Depending on deformation in a cell the calculated jump in displacement can reflect an opening mode (u n > 0) or a closing mode (u n < 0). An opening mode is initialized as a crack with the amount of opening prescribed by the kinematic analysis.
The discontinuities in a closing mode are initialized as cracks with the calculated direction and tangential jump in displacement, but with zero normal jump in displacement. After the initialization, the simulation is run for 15 more days. The initial fractures may open or close, and other fractures may be created and evolve. The resulting predicted values of u n , u t and α are compared to the kinematic analysis of RGPS data.

Statistical Method
We propose a novel emulation and calibration framework based on proper distance metrics defined on the image space, this allows us to properly account for the geometric differences between model outputs and observations. In this section, we first introduce the concept of phase and amplitude distances provided by image warping, where the first one quantifies the amount of geometric deformation required to correct for misaligned and misshapen model errors, while the latter quantifies the remaining model error after optimal alignment. We then describe how to construct an emulator and carry out calibration based on these two scalar metrics. This simplifies the high-dimensional image model outputs to scalar responses and significantly reduces the complexity of the statistical emulator and parameter inference procedure.

Image Warping
Our approach to image warping is to use the concept of phase and amplitude of an image that was introduced by [Xie and Srivastava, 2016]. This work has been motivated by the the same concept developed in the work on shape analysis of objects and phase-amplitude modes in functional data (see Srivastava and Klassen [2016], Marron et al. [2015], and Tucker et al. [2013] and references therein for a review of the concept). A deformation is defined to be a change of only the phase of an image (where the pixels are located) and not the amplitude. This view point then sets up the registration problem.
Let the set of all images be F = {f : D → R n |f ∈ C ∞ (D)} of the images in the domain D that map to pixel values in R n . Here we require n ≥ dim(D) for the theory presented in this section. An example of a suitable image is that of an RGB image. In this case, D = R 2 and n = 3, the space of red, green, and black intensities at a given s 1 , s 2 ∈ D. In the case of a gray scale image, where f (s 1 , s 2 ) : D ⊂ R 2 → R, one may instead consider an analogous image using the gradient of f such that F (s 1 , s 2 ) = (∂ s 1 f, ∂ s 2 f ) or a combination of intensity and gradient such as F (s 1 , s 2 ) = (f, ∂ s 1 f, ∂ s 2 f ). We discuss our particular choices for image type in Section 4.2. Also let Γ = Diff + (D) be the set of of all orientation-preserving diffeomorphisms that preserve the boundary of D. It can be shown that Γ is a group that contains an identity and inverse element. The registration of two images f 1 and f 2 is expressed as finding γ * such that (1) is not a proper Riemannian metric; specifically, it lacks the symmetry property.
To overcome this problem Xie and Srivastava [2016] introduced the q-map, where for an f ∈ F, we have a mapping Q : F → L 2 for any s ∈ D, where a(s) = ||Jf (s)|| V and Jf (s) is the Jacobian matrix of f at s. The term a(s) is known as the generalized area multiplication factor of f at s. The use of the q-map provides a distance that is a proper metric on the space of images for the registration problem For proofs of the metric properties of (3) the reader is referred to Xie and Srivastava [2016]. This registration problem can then be solved using gradient descent.
The gradient descent is computed iteratively over Γ, where at each iteration the directional derivative is taken on the tangent space of Γ at the identity element (T γ id (Γ)). The directional derivative is computed using an orthonormal basis expansion on T γ id (Γ) and details are given in Xie and Srivastava [2016]. The current implementation of the gradient decent scheme is sensitive to step size and the number of basis functions used to represent γ * in the orthornomal expansion. These can be tuned, however, for the specific problem.
This registration problem then gives rise to the definition of amplitude distance between two images We can then also define a distance on the amount of deformation required to register the images.
Since the space Γ does not have a simple Riemannian geometry or one cannot be found through a transformation as in [Tucker et al., 2013] we will use an extrinsic approach. We will use the q-map of the warping function and the resulting L 2 metric. This gives us a nice metric for computing the distance, but further statistical analysis becomes hard, such as principal component analysis, due to the fact that the inverse mapping is not available in closed form.
We will call this the phase distance and it is defined as where q γ = Q(γ). This distance measures how far the computed warping is from the identity warping. We will use these two distances in the calibration solution discussed later on.

Statistical Emulator
The sea ice model represents the formation and evolution of leads(cracks) resulting from the complex interactions among input processes. Due to the computational cost, we have a limited number of model runs. A standard approach is to replace the computer model by a fast surrogate statistical model often called an emulator. Typically, the first step is to obtain model outputs at a set of design points, then construct a statistical model such as a Gaussian process to interpolate these model outputs [Oakley and O'Hagan, 2002, Higdon et al., 2004, Conti et al., 2009.
Instead of interpolating model outputs, we use an emulator to interpolate the amplitude and phase distances, d a and d p , provided by image warping. We let Y(θ) ∈ F denote the model output corresponding to input parameter θ ∈ Θ ⊂ R d for an integer d ≥ 1, then the dis- We construct a Gaussian process emulator η(θ, d) with mean µ = (µ a , µ p ) T and covariance function C(·, ·) for the distance metrics d(θ). For the mean we use a linear function µ k = x T β k for both k = a and p, where x are the covariates and defined as the input parameter θ plus an intercept term, and β k is a d + 1-dimensional vector of regression coefficients. The covariance function models the dependence between every pair of metrics in terms of their distance in the parameter space. For example, we have used the exponential covariance function in our simulated example and data application. It is a commonly-used, separable, and stationary covariance function with the following form, where ρ k1 , . . . , ρ kd are the range parameters, σ 2 k is the partial sill, and τ 2 k is the nugget. The range parameters control the rate of the correlation decay with distance for each dimension of the input space. Here we use ξ k = (ρ k , σ 2 k , τ 2 k ) to denote the parameters for the statistical emulator, to distinguish them from the physical model parameters θ. We consider independent Gaussian processes for amplitude and phase distance since these two metrics are uncorrelated and hence we set Cov (d a (θ), d p (θ )) = 0. Other covariance funtions can be used here; we suggest performing cross-validation to select the appropriate covariance function.
The emulator is then constructed based on the metrics evaluated at the chosen design points using Latin hypercube sampling (LHS) [McKay et al., 1979]. Let X be an N × (d + 1) covariate matrix and d k = (d k (θ 1 ), . . . , d k (θ N )) T , k = a, p be the vectors of amplitude and phase distance metrics evaluated at the design points. The Gaussian process emulator described above specifies d a and d p as the following independent multivariate normal random variables: We estimate the emulator parameters β k and ξ k by maximizing the likelihood given by the above model and denote the resulting maximum likelihood estimators (MLE's) asβ k andξ k . For any new input parameter setting θ ∈ Θ, the emulator predicts its corresponding distance metrics through the conditional normal distribution of d a (θ) | d a and d p (θ) | d p , using the MLE'ŝ β k andξ k as plug-in estimates. This predictive process serves as the surrogate model in the following section on model calibration.

Model Calibration
We use Z d = (d a , d p ) T = (0, 0) T as the observation for calibration, since we use the observed image Z as the reference, and registering an image to itself results in zero amplitude distance and zero phase distance. We assume that the data model for Z d is where θ * is the 'best' parameter setting that yields d(θ * ) = 0. The data-model discrepancy = ( a , p ) T is assumed to follow a truncated normal distribution with the following density: with the discrepancy standard deviation ψ k > 0 for k = a, p and the probability and cumulative density functions φ and Φ for the standard normal distribution. We have chosen a truncated normal distribution with lower bound zero since we assume that the distance metrics for model outputs cannot take lower values than those for the observational data Z.
We choose to adapt the standard Bayesian calibration framework for parameter inference as it provides more detailed information on input parameter uncertainty, in the form of posterior densities [Kennedy and O'Hagan, 2001]. Each input parameter in θ receives a uniform prior with lower and upper bounds determined by its physical limits. Each discrepancy variance parameter ψ 2 k receives an informative inverse gamma prior with a shape parameter a k = 20 and a scale parameter b k = a k d 2 k,10 th for k = a, p, where d k,10 th denotes the 10 th quantile of the distance metrics d k obtained at the design points. This prior encourages the resulting full conditional posterior for the standard deviation ψ k to have its mode around d k,10th , so that roughly 10 percent of the existing model runs at the design points are covered by the 95% credible interval. Informative priors are chosen for the variance parameter because we have only one observation for calibration; this inverse problem is non-identifiable without strong priors. The posterior samples for θ, ψ 2 a and ψ 2 p are obtained from Markov chain Monte Carlo using the Metropolis-Hastings algorithm.

Experiments
We conduct two experiments to illustrate the usefulness of the proposed method. In the first experiment, we generate model outputs from applying a diffeomorphism on a template that represents a feature similar to a crack. The second experiment is based on the sea ice model outputs, where one run of the model is randomly chosen to be the synthetic observation. The experiments establish proof of concept, the first in a simplified setting and the second with a complex geometry that is consistent with the real data application. The experiments also establish the applicability of this method to discrete data sets, even though the theory itself is developed in a space of smooth functions. In both experiments, we compare calibration results using the Euclidean distance, phase distance, amplitude distance, and both phase and amplitude distance combined.

Diffeomorphism Warp
We illustrate the calibration method using a toy example with geometric features similar to the ice opening features resulting from ice deformation. The mathematical model used here applies a boundary preserving differmorphism to a 2-dimensional surface on the unit square, with input parameters θ = (θ 1 , θ 2 ) to define the map. This allows us to evaluate our calibration method on a data set where simulation differs from observation by an actual diffeomorphism, an ideal case. Let s = (s 1 , s 2 ) denote the coordinate of an image pixel. The diffeomorphism W : R 2 → R 2 maps a 2-dimensional coordinate to a new coordinate that is also 2-dimensional.

(9)
We generated 50 model runs at a set of design points chosen using Latin hypercube sampling. Figure 1a shows the template image for generating the model runs. For each design point, the template is deformed to create a warped image based on equation (9). As an illustration, Figure 1 shows the model outputs for four different input parameter settings. The branch located at the bottom right corner of the template image is stretched to the right when θ 1 is positive and compressed to the left when negative, while a similar deformation action is applied to the vertical direction controlled by θ 2 . We also generate a synthetic observation (Figure 1b) from the diffeomorphism warp using input parameter (0.3, 0.1). Therefore, all the model runs here have misaligned and misshapen features compared to the synthetic observation.
We compute the deformation metrics for each model run using the synthetic observation as the reference. In this case our image is of the form F (s 1 , s 2 s 2 ) is the pixel intensity. We calculate the gradient using finite difference approximations.
The optimal deformations for the four example runs and the warped model runs are shown in Figure 1. We see that all four warped model runs are indistinguishable to the synthetic observation visually. This indicates the deformation has been applied correctly and estimated well.
We fit independent GP emulators to each metric then perform Bayesian calibration as described in Section 3.3. In our MCMC implementation, we sample three 30,000 posterior samples at different initial values with 15,000 burn-in; MCMC convergence is assessed using the Gelman-Rubin diagnostic [Plummer et al., 2006]. Calibration results using different metrics are presented in Figure 2. It is notable that the posterior distributions from calibration using the amplitude or Euclidean alone are very flat, providing little information on the input parameter, while the posterior using the phase distance peaks near the true input parameter value, and the calibration using the combined metrics is similar to using phase dis- The results confirm that the proposed method is able to properly account for misaligned and misshapen model outputs. Because the model runs and the synthetic observation for this toy example are all warped from the same template, all model runs have the same amplitude distance except for small numerical differences due to interpolation in generating the images.
This means there exists a γ * that exactly maps each image to the synthetic observation, therefore, the phase distance quantifies the amount of warping in order to correct for misaligned and misshapen features. The calibration using the phase distance alone suggests the input parameter setting that produces the correct geometric features which are least misaligned and misshapen. This is particularly useful for model development and evaluation as a measure of a model's ability to produce correct coherent structures, a geometric feature that persists in a flow field. Moreover, for models which produce highly localized coherent structures, hurricanes for instance, the Euclidean metric may completely disqualify a model run which produced the correct hurricane characteristics just in the wrong location. In this case the model itself may be well calibrated with location error coming from external sources, such as error in an initial condition taken from observations.

Sea Ice Model Experiment
In the previous toy example our model runs were generated by an actual diffeomorphism, however in most geophysical applications model output comes from a physical numerical model. These model outputs corresponding to different parameters may not be an exact diffeomorphism warp of each other. Here we perform a perfect model experiment for the sea ice model described in Section 2. In this perfect model experiment we take the observation to be the model output for a chosen set of parameters. This experiment demonstrates, through the recovery of the chosen model parameterization using the calibration scheme, that the method works well for physical models whose output is not generated from a diffeomorphism.
We produce 80 models runs for a set of To carry out the analysis, we compute the optimal warping function for each model run.
Assessing the metrics corresponding to the model runs, we notice that the phase distances are left skewed. As a result, we compare cross-validation performance on GP emulators fitted to the original and the log transformed distances; results indicate that data transformation is necessary for this experiment to ensure the surrogate model has good prediction performance, for instance, small mean square prediction error and reasonable coverage.
Bayesian calibration and the necessary convergence diagnostic are then carried out; the posterior density functions are shown in Figure 3. The marginal posterior mode and credible interval for different metrics are summarized in Table 1

Application to Sea Ice Model and Sattelite Data
The proposed method is applied to the RGPS data discussed in Section 2. We warp each model output to RGPS data using the same choice for step size and the number of basis functions to represent γ as in the sea ice experiment in Section 4.2. MCMC convergence has been assessed using the Gelman -Rubin diagnostic [Plummer et al., 2006]  and discarding the first 15,000 for burn-in. Calibration results are presented in Figure 4. The posterior modes for the tensile and shear strength parameters from using the combined metrics are 25.45 (11.601, 48.063) and 75.71 (60.514, 95.42), respectively, with the credible intervals given in parenthesis. To establish that the posterior mode represents a viable parameterization, we run the elastic-decohesive model and assess the output with input parameters equal to the posterior mode, shown in Figure 5. In contrast to the perfect model case above, we cannot know whether there is a set of parameters which can exactly reproduce the RGPS data, in fact likely not. To quantify the performance of the model with the calibrated parameters we compare both the Euclidean and q-map distances (i.e. L 2 distance between q-maps of the images) between the model runs and the RGPS data, for both the calibrated parameter settings and the original parameter settings chosen from physical considerations in . We also apply our warping algorithm using each model output to measure the amplitude and phase distances between the runs and the RGPS data. Before preforming any warping we find that the parameter settings used in  give q-map and Euclidean distances of 23.9 and 46.3, respectively, between the model run and observations compared to that of 18.2 and 42.4 for the model run using the calibrated parameters. The amplitude distances for both model runs are similar at 11.2 for the original parameters and 11.9 for the calibrated parameters.
The most striking difference comes from the phase distances at 5.4 for the original parameters and 1.1 for the calibrated parameters. This suggests much less warping was required to align geometric features when comparing the calibrated model run to the observations. We believe the calibration based on both amplitude and phase provide a more reliable inference on input parameters than using Euclidean alone based on the perfect model study where the truth is known. In the instance shown in Figure 3 the Euclidean distance is less informative with a very diffuse posterior and a mode which misses the true value. Further, the center mode in Figure 4 encompasses that chosen from physical considerations in Sulsky . In the Appendix we show the case where our image was taken to be f (s) = (|[u]|(s), α(s)).

Discussion
The calibration method presented is the first of its kind based on the image warping metrics.  is of importance. Further, the metric itself is a proper distance provided the optimal γ is found.
Currently, we are using gradient descent, local minima are a possibility. The gradient decent scheme can be very sensitive to step size and the number of basis functions chosen to represent γ. In this work, we manually tuned the step size and number of basis functions to produce an optimal result. Future work would include looking at other optimization methods. The proposed method also reduces the high-dimensional model outputs into scalar responses, such that We illustrated the usefulness of the proposed method through a toy example and a perfect model experiment using the elastic-decohesive model. We used 50 design points for the first experiment and 80 for the latter. While the number of design points required depends on the complexity of the response surface and the number of active inputs influencing the response, the sample size we used here was sufficient for our experiments. A rule of thumb for selecting the number of design points is 10d, where d is the dimension of the input space [Loeppky et al., 2009, Harari et al., 2018. While the results in our experiments show that the proposed method performs well, we note that we have not added observation error when generating the synthetic observation. This can be a major caveat where a more rigorous investigation should be conducted. It remain future work to study how observation error could affect warping and calibration results.
The ultimate application of this method should be to real world observations. The use of image warping techniques is motivated by a need to assess sea ice model output consisting of fractures. For climate modeling, it is not necessary to get these features exactly right. In fact, it is not likely that we would be able to predict the cracks exactly since the forcing that drives the motion in the model is too uncertain. Thus, we desire a metric that can compare data and model favorably if the predicted cracks can be warped to the observed cracks that is, if the predicted cracks are somewhat misaligned or misshapen. In our study using the satellitederived RGPS data of sea ice deformation we are able to obtain informative posteriors using the combination of amplitude and phase metrics, and the posterior mode seems to represents a viable parameterization, when compared to that chosen from physical considerations in Sulsky   Magenta color "+" indicates the input parameters corresponding to the synthetic observation.
Bottom: RGPS data application: Posterior density function of input parameters using both amplitude and phase, amplitude only, and phase only.