Assessing ice sheet models against the landform record: The Likelihood of Accordant Lineations Analysis (LALA) tool

Palaeo‐ice sheets leave behind a rich database regarding their past behaviour, recorded in the landscape in the form of glacial geomorphology. The most numerous landform created by these ice sheets are subglacial lineations, which generate snapshots of the direction of ice flow at fixed (yet typically unknown) points in time. Despite their relative density within the landform record, the information provided by subglacial lineations is currently underutilised in tests of numerical ice sheet models. To some extent, this is a consequence of ongoing debate regarding lineation formation, but predominantly, it reflects the lack of rigorous model‐data comparison techniques that would enable lineation information to be properly integrated. Here, we present the Likelihood of Accordant Lineations Analysis (LALA) tool. LALA provides a statistically rigorous measure of the log‐likelihood of a supplied ice sheet simulation through comparison of simulation output with both the location and direction of observed lineations. Given an ensemble of ice sheet simulations, LALA provides a formal, and statistically underpinned, quantitative assessment of each simulation's quality‐of‐fit to mapped lineations. This enables a comparison of each simulation's relative plausibility, including identification of the most likely ice sheet simulations amongst the ensemble. This is achieved by modelling lineation formation as a marked Poisson point process and comparison of observed to modelled flow directions using the von Mises distribution. LALA is flexible—users can adapt parameters to account for differing assumptions regarding lineation formation, and for variations in the level of precision required for differing model‐data comparison experiments. We provide guidelines and rationale for assigning parameter values, including an assessment of the variability between users when mapping lineations. Finally, we demonstrate the utility of LALA through application to an ensemble of simulations of the last British‐Irish Ice Sheet. This comparison highlights the benefits of LALA over previous tools and demonstrates some of the considerations of experimental design required when identifying the fit between ice sheet model simulations and the landform record.

The information left behind by palaeo-ice sheets provides a long-term record of ice sheet behaviour (Clark et al., 2018), which, if unlocked, has great potential to be used to improve the predictive ice sheet models essential to forecasting the fate of ice sheets in our warming world (e.g., Edwards et al., 2021;Goelzer et al., 2013;Lipscomb et al., 2021;Nowicki et al., 2020), and our reconstructions of past ice sheets (e.g., Andrews, 1982;Stokes et al., 2015).Quantitative comparisons between numerical model simulations and palaeo-ice sheet evidence are becoming more commonplace (e.g., Ely et al., 2021;Gandy et al., 2019Gandy et al., , 2021;;Tarasov et al., 2012), replacing qualitative descriptions of fit to evidence (e.g., Boulton & Hagdorn, 2006;Siegert et al., 2001).This shift is prompted by a need to quantify the degree-of-fit between models and evidence to assess output uncertainty and is facilitated by an increase in computing power enabling sufficient resolution within simulations for comparison to take place.However, the use of quantitative model-data comparison tools in palaeo-ice sheet modelling experiments are far from routine.This is at least partially a consequence of the underdevelopment of model-data comparison tools.
The most abundant landform evidence left behind by palaeo-ice sheets are subglacial lineations (henceforth lineations).Often thought of as the separate categories of drumlins (e.g., Menzies, 1984), crag and tails (e.g., Dowdeswell et al., 2016), and mega-scale glacial lineations (MSGL) (Clark, 1993), lineations can be broadly defined as streamlined hills, formed at the ice-bed interface.Lineations are typically on the order of 100-1000 m in length, though they can reach several kilometres in the case of MSGL (Spagnolo et al., 2014).The origin of subglacial lineations, especially drumlins, has garnered much scientific debate and is an active field of research (see, e.g., Ely et al., 2022).However, general agreement can be found amongst geomorphologists on two factors: (i) Lineations are streamlined in the direction of ice flow, and thus represent former ice flow direction at a point in time during the lifecycle of an ice sheet; and (ii) lineations form under warm-based conditions, conducive to the transport of subglacial material.Two further observations help glacial geomorphologists reconstruct the past behaviour of ice sheets.First, lineations typically occur in fields of regular arrangement with similar orientation and morphology (Clark et al., 2018), although isolated examples do also exist (Evans et al., 2015).As such, lineations can be grouped into flowsets, larger regions of a palaeo-ice sheet bed thought to represent an ice sheet flow event (Clark, 1997).Second, sets of lineations can be observed to cross-cut each other, with younger forms superimposed upon those formed during older flow events (Clark, 1993).
Such cross-cutting relationships enable a sequence of ice flow orientations to be deciphered.On this basis, the interpretation of lineations has provided much insight into the operation of palaeo-ice sheets (e.g., Dyke, 2008;Greenwood & Clark, 2009a;Hughes et al., 2014;Stokes et al.2009).
Despite being a pervasive and information-rich source of data on palaeo-ice sheet behaviour, lineations are underutilised in palaeo-ice sheet modelling experiments.To date, we are not aware of any statistically rigorous approach for incorporating observed lineations into model-data comparisons.A first attempt at building a tool for comparing simulated ice flow directions and those derived from observations of lineations has however been developed by Li et al. (2007).The Automated Flow Direction Analysis (AFDA) tool provides a measure of the degree-of-fit between simulated and observed flow directions throughout the duration of an ice sheet simulation (Li et al., 2007).
Two metrics are calculated from gridded datasets of simulated and observed flow directions.The resultant mean difference aims to provide a measure of the overall directional offset between simulation and data, whilst residual variance is used as a measure of the level of agreement between the shape of the two flows.Ely et al. (2021) subsequently developed a workflow for assessing whether cross-cutting relationships were replicated, whereby model-data agreement was declared if a simulation was able to replicate flow directions in the correct sequence.
Despite its existence, uptake of AFDA has been low.As of February 2023, a Google Scholar search indicated that Li et al. (2007) had been actively used in six studies that use the tool within ice sheet model experiments.This may reflect the tool being ahead of its timeice sheet simulations have only recently been able to simulate the detailed flow fields recorded in lineations and written in the now largely defunct, GIS-based language.It may also be due to certain weaknesses of AFDA, which we aim to address in this work.The lack of a comprehensive and formal statistical underpinning to AFDA makes decisions about the degree-of-fit between an ice sheet simulation and observed lineation data highly subjective-the declaration of a model-data match requires user-defined thresholds in resultant mean difference and residual variance (Ely et al., 2021).This hinders comparison between simulations.Furthermore, there is a growing realisation within ice sheet modelling that emulators, statistical surrogates of numerical models, may be useful for cutting computational cost and exploring unsampled parameter space (e.g., Edwards et al., 2021;Tarasov et al., 2012).It may be desirable in the future to emulate the output of a model-data comparison tool, to identify the parameter-space most likely to replicate flow directions recorded in lineations.The binary fit/no-fit metric of AFDA is unsuitable for emulation which requires a more refined, continuous, and rigorous measure of fit.
In this paper, we present a new, statistically rigorous, tool for comparison of ice sheet simulations with observed flowsets.Our Likelihood of Accordant Lineations Analysis (LALA) tool provides a quantitative, and continuous, measure of the likelihood of the lineations given an ice sheet simulation.LALA has a formal, and robust, methodological underpinning.We develop a probabilistic model linking an ice sheet simulation to the creation of lineations using a marked Poisson point process (Kingman, 1993).In addition to comparing observed flow direction with the ice sheet simulation, LALA also incorporates the number and locations of the observed flowsets.In doing so, key information regarding both the presence and absence of lineation evidence in certain locations is incorporated into the model-data comparison.The LALA tool is designed to provide an objective assessment of model-data fit yet is also highly adaptable to allow incorporation of expert knowledge regarding lineation formation.
Our paper is laid out as follows.In Section 3, we provide an intuitive explanation of the underpinning probabilistic model that we use to link an ice sheet simulation to observed flowsets.This probabilistic model is constructed from three simple assumptions.We also describe how our LALA tool follows directly from these initial assumptions.After this key mathematical framework is demonstrated, a worked toy example is presented in Section 3 to highlight the premise behind the tool in a less-mathematically taxing manner.A critical parameter when considering the model-data comparison is the uncertainty in the ice flow orientation derived from interpreting lineations.Section 3 presents results from an experiment designed to quantify this uncertainty.Finally, in Section 3, we demonstrate the utility of the tool by applying it to a small ensemble of simulations of the last British-Irish Ice Sheet (BIIS).

| CONSTRUCTION OF THE LALA TOOL
Our tool, LALA, is provided as a Python script (v3.0+) along with a tutorial that works through synthetic examples of observed lineation locations and directions, and ice sheet model output (see Appendix A and Data Availability Statement).It is built to handle netCDF files (a common format for ice sheet model output) of ice sheet model simulations and observations of lineation formation.An overview of the statistical underpinning of LALA is provided below, with a fuller description provided in Appendix A. In brief, LALA considers that lineations formed under specific subglacial conditions, and aligned with the ice-flow direction at the time of their formation.The likelihood that the observed lineations were formed by a prescribed ice sheet simulation is calculated, providing an assessment of that specific simulation's quality-of-fit.Typically, we expect our model will be used to compare different simulations within a large ensemble, perhaps representing a range of hypothetical forcings and ice-sheet parameter selections, to identify those simulations that are most likely and hence narrow down the plausible range of forcings and parameters.The output from LALA can however also be considered on a per-flowset basis, as demonstrated later (Section 3).The relative timings of lineations can be inferred from the observational record and are a valuable addition to learning about the timings of past flow directions of palaeo-ice sheets.Accounting for these so-called cross-cuttings is currently beyond the scope of LALA, but would be an interesting extension to explore.In this paper, we only consider a point in the centre of each flowset, which may not represent flow across complex flowsets with curving directions.Equally, users of LALA could score each model cell a flowset covers and check for coherent temporal matches.

| Scoring concept-A simple probabilistic model for flowset formation likelihood
LALA is built upon a rigorous statistical foundation: A simple probabilistic model, based upon three initial assumptions regarding the formation of lineations, that provides a link between an ice sheet and its resultant lineations.Having constructed this probabilistic model, the LALA tool and the score for any ice sheet simulation follows directly.
We suppose that we have observations of lineations, gridded to the scale of the ice sheet model output we wish to compare, which correspond to n flowsets across our overall study region X, a complete or partial palaeo-ice sheet bed.These flowsets are observed at locations x 1 , …, x n and are accompanied by estimates of their inferred directions θ 1 ,…, θ n .Given these observed flowsets, providing paired ðx i , θ i Þ n i¼1 information, we aim to score any hypothetical/simulated ice sheet M to assess its level of agreement with the observations.This score is obtained by evaluating the log-likelihood of the ice sheet under con- We will denote an ice sheet at time t and location x as Mðx, tÞ.
This Mðx, tÞ might consist of many variables, for example, ice thickness, basal and englacial velocity (speed and direction), thermal regime, and mass-balance.Given an ice sheet Mðx, tÞ, over locations x X and times t T , we assume flowsets are created according to the following three basic principles: 1. Ice cover is required for lineation formation.At any time t when a location x is covered by ice, there is a small probability that lineations could form in that location.This probability of forming lineations may depend upon the location, the time, and the properties of the ice sheet at that time.We denote the probability of lineation formation at location x and time t as λðx,t, Mðx, tÞÞ.This will critically depend upon the properties of the ice sheet under consideration and the location x.
2. The formation of lineations are completely independent from one another in both space and time (conditional on the value of λ).This assumption likely holds at the scale of an ice sheet that we are concerned with here, though we note that at the scale of neighbouring lineations, interactions between lineations may occur (Ely et al., 2018).
3. Lineations align with ice flow direction.Intuitively, the lineations will align (at least approximately) with the ice flow at the point they are formed.However, from lineation morphology alone, it can be difficult to ascertain the upstream and downstream ends of a lineation (Spagnolo et al., 2010(Spagnolo et al., , 2011)), especially in regions that experienced a complex ice flow history.Thus, we allow for orientations to be exactly opposite of those prescribed by the user (i.e., the lineations could record ice flow in the exact opposite direction).This assumption could be relaxed in future versions of LALA, to account for regions where ice flow direction is well known.
Under these initial assumptions, the flowsets form what is known as a marked, inhomogeneous, Poisson point process (Kingman, 1993)-with the mark being the direction of the flowsets.A Poisson point process is used to model points, considered to be independent, randomly over a set space with a specified rate.A marked Poisson point process is the same idea, with the added benefit of adding additional information to each point, in this case: direction.
Here, lineations and associated orientations are the points, and the model domain is the space we are considering.The associated rate that is attached to this Poisson process is calculated from the expected number of lineations in the domain over time divided by the number of time-integrated plausible areas for lineation formation, and the inhomogeneous property allows for the rate to vary with respect to space and time.
Given our set of n observed flowsets at locations x 1 , …, x n and with estimated paired directions θ 1 , …, θ n , we can use this Poisson point process model to calculate the log-likelihood of any particular ice sheet.This log-likelihood is the score of our LALA tool.Those ice sheets with greater log-likelihoods are more likely than those ice sheets with lower log-likelihoods in the sense that they make the observed lineations (in terms of both location and direction) more probable.The full details on constructing the log-likelihood can be found in the Appendix but we provide an intuitive explanation of the main elements below.
We restrict our study region to those areas which have been assessed for lineations (i.e., we are aware whether the location has, or does not have, lineations).Those regions that have not been mapped (and hence where the existence, or absence, of lineations is unknown at an ice sheet model resolution) should not be included in our tool or form part of the study area X .When scoring an ice sheet against the flowsets, the value of the ice sheet in this unmapped area should also be discarded.

| Overview of scoring direction and location
The rate of lineation formation λðx, t, Mðx,tÞÞ can, in principle, depend upon multiple variables: The location, the time, and the properties of the ice sheet and properties of the sediment.However, for initial implementation of our tool and to simplify the intuitive explanation below, we will reduce the dependence on these other variables so that λðx, t, Mðx, tÞÞ has two states: (i) Potential for lineation formation and (ii) impossible for lineation formation.We set λðx, t,Mðx, tÞÞ ¼ 0 when the location x is not grounded ice, the ice thickness is less than 10 m, or the speed of the ice sheet flow is below 10 ms À1 .At such times or locations, lineation formation is considered impossible.We focus on these simplistic rules here; however, these assumptions can be updated to account for other information considered by the user to be important for lineation formation if available from the ice sheet simulation.Otherwise, when the ice sheet has lineation formation potential we will choose a constant λðx, t, Mðx,tÞÞ ¼ λ.This assumes that where and when lineation formation can occur, it will do so at a similar rate across the ice sheet bed.Practically, this means each time step will have a map of where lineations can and cannot form, according to the above criteria.Where formation is possible, we assign a rate λ and where formation is impossible assign a rate 0. The more flexible version (with a general λ) is described in Appendix A.

| Number and locations of flowsets
The first two of our assumptions (ice cover and the complete independence of lineation forming events) imply that for any ice sheet M, the number and locations of the flowsets follow an inhomogeneous Poisson point process (Kingman, 1993).We can calculate the likelihood of not only observing n flowsets but also that they fall in the locations x 1 , …, x n where they have been observed: Here, in our simplified two-state λ case, the term Λ M ðx i Þ reduces to λT M ðx i Þ, where T M ðx i Þ is the total length of time that location x i has the potential to form flowsets, that is, is suitably covered by the ice sheet M. The term Ð X Λ M ðx i Þdx ¼ Λ M ðX, TÞ integrates these individual values over the whole study area and effectively calculates the total "time Â area" when the ice sheet M has the potential to form lineations.
Note that under our model, the number of flowsets for a particular ice sheet M follows a Poisson distribution with mean Given the same lineation formation rate λ, larger ice sheets (covering greater areas where lineation formation is possible) will be expected to create more lineations than smaller ice sheets.Additionally, the inclusion of the term 1) highlights that we are more likely to see lineation formation in locations where the ice sheet has remained for longer.This will be reflected in the scores that we give to different ice sheet simulations.Solely in terms of scoring the location of the lineations, we will tend to give greater scores to those models where the ice sheet is concentrated on the lineation locations; and penalises simulations that extend beyond the assumed glaciated area through time ðΛ M ðX ,TÞÞ, as defined by the user.

| Flowset orientation
In addition to considering the fit to the locations of the observed lineations, we also seek to assess the fit between the orientations of those lineations and the simulated ice sheet.Suppose that there is a flowset which forms at location x i and time t ?i .We assume that the orientation of the resultant flowset aligns with the direction of ice flow in the location at the point of formation.Specifically, we model the observed orientations θ i according to a von Mises distribution (von Mises, 1981) centred around μðx i , t ?i Þ, the ice flow direction at the (unknown) time of formation for flowset i.This von Mises distribution provides a distribution for angles that is an approximate analogue of the normal distribution.Due to the morphology of lineations, we account for inferred observed directions occurring in the reported direction or the exact opposite direction.This leads to us modelling the observed direction θ i as where In reality, we do not know the precise time t ?
i at which the lineations were formed in location x i , only that it must have occurred at one of the times when the ice sheet had lineation formation potential.
To assess the overall directional fit, we must therefore average over all these times when lineations could be formed.Precise details are given in Appendix A, but, in our simplified two-state case for λ, the likelihood of a particular ice sheet model M in terms of the agreement between the observed direction of lineation i and the ice sheet is therefore: Here, μðx; tÞ is the orientation of ice flow in model M in location x at time t and f M ðθ i ; μ, κÞ is the mixture of the von Mises and flipped von Mises with location μ and concentration κ as described in Equation ( 2).
We combine these directional components, multiplying the likelihood for each individual flowset, over all the n flowsets within LALA.
Intuitively, this directional likelihood component to LALA will provide greater scores to those ice sheets where, when lineation formation is possible at location x i , the ice sheet flow aligns precisely with the observed direction of the lineations.Low scores will be obtained for those ice sheets that, when lineation formation is possible, have flow directions that lie at odds with the observed direction of the flowsets.

| Final ice sheet log-likelihood lðMjfx
To calculate the log-likelihood of any hypothetical ice sheet M given the number, locations, and directions of the flowsets observed in the study region X , we combine the components described above and take logs.In the case of our simplified two-state lineation formation model, and excluding the possibility that a lineation may have formed outside the period of study, the final log-likelihood for a simulation M becomes The first component of this log-likelihood (within the square brackets of Equation 3), relating to the fit of the simulated ice flow direction compared to the flowset orientation, can be used as an intuitive indicator for which flowsets are contributing the most to the final score.We can use this as a way to judge the flowsets, which fit well to a simulation and which do not.This idea is explored further in Section 3.

| Accounting for lineations that occurred outside the simulation period
We anticipate that LALA may be used for a wide range of modeldata comparison experiments, with different time periods and areas.
As such, some study areas may contain glacial lineations that were formed outside of the time period being simulated by the ice sheet model (e.g., during a prior glacial or later/earlier during the same glaciation).Additionally, it is possible that lineations may have formed as a result of processes that are poorly represented in the ice sheet model, or there may be insufficiently high temporal resolution of the simulation output to capture all lineation forming events-that is, lineation formation occurred between simulation time steps.
Without additional interpretation of the evidence, we cannot rule out such possibilities.This possibility is accounted for in LALA by assigning a small, but statistically relevant, probability that any observed lineation may be unrelated to the ice sheet simulation we wish to test.
We can incorporate all of the above possibilities into our loglikelihood tool by adding an additional component to our Poisson point process.This additional component assumes that, at some time before the ice sheet simulation we wish to test starts, any location x may have already had a lineation forming event.If such a prestudy event has occurred, then we presume that the resultant lineation may be aligned in any direction uniformly.Users are able to specify the probability, p, that any given lineation relates to such a pre-study ice sheet (i.e., it is unrelated to the simulation period we are testing).To account for this in the scoring, we extend our Poisson process model to incorporate a prestudy pseudo-time period during which we assume a different rate of lineation-forming λ ?(chosen adaptively according to our selected probability p).Lineations created in this pseudo-time pre-study period are modelled to align uniformly in any direction.Further details are given in Appendix A and an explanation of how to find the appropriate value λ ? is given below.

| Parameter selection
The use of LALA requires several tool parameters to be determined.
Our tool allows users to select suitable values based upon their expert knowledge.However, we also propose automated choices as follows: λðx,t, Mðx, tÞÞ is the rate of lineation formation during the study period.The choice of λðx, t,Mðx, tÞÞ is critical to LALA's final loglikelihood.In the absence of detailed expert information, we suggest the two-state approach described above whereby we assume lineation formation is impossible (i.e., λðx, t, Mðx, tÞÞ ¼ 0) if the location x is not grounded; if the ice sheet does not cover location x at time t; or if the ice sheet has a flow speed below 10 ms À1 or a thickness of less than 10 m.Otherwise, we assume lineation formation is possible and λðx, t, Mðx, tÞÞ ¼ λ.The conditions for the inability of lineation expected number will depend upon λ.We then select λ so that E M † ½N equals n, the number of flowsets actually observed in our region of study.This plausible value λ is then fixed and applied for the scoring of all ice sheets simulations under consideration by the tool.
λ ?ðxÞ is the rate of lineation formation outside the study period.
We also consider in a similar way a rate for lineation formation occurring before the study period.We create a grid matching the simulation domain area that indicates areas where conditions were conducive to lineation formation.For example, users may want to consider limiting to areas of the palaeo-ice sheet bed where there is adequate sediment available for lineation formation.Furthermore, regions of the simulation domain that cover the deep ocean cannot host lineations as an ice sheet is unable to ground there.To rule out these areas where lineation formation is impossible, we set λ ?ðxÞ ¼ 0 where formation is impossible and λ ?ðxÞ ¼ λ ?elsewhere.This rate also depends on the user-defined probability p. p is the probability that an observed lineation relates to a period of time outside the simulation period.The user is able to specify the value of p as they see fit; however, the default for LALA is to select p ¼ 0:01 (i.e., out of 100 observed flowsets, we might expect 1 to have formed outside of the realm of the studied time period).
Note that all of the parameters above are dimensionless.When using this tool, the values of κ and p can be changed according however the user wants to use the tool.

| TOY EXAMPLE
To demonstrate the principles of LALA, we here present a toy example.We consider a simulation, say M, with a 5 grid cells by 5 grid cells, one flowset and three time steps.In Figure 2 actually conducted in radians (Section 3).LALA will perform this transformation from degrees to radians for the user if the input variable for radians is set to False.For this example, we use κ ¼ 5.
First, we need to calculate values of λ and λ ?, the rate of lineation formation during the study and in our prestudy, respectively, using a simulation within the ensemble.This simulation, say M † , can be identified through comparison to other metrics, such as ice extent (see Section 3).We assume that in the prestudy any grid cell could have have expected where E M † ½N is the expected number of flowsets formed for simulation M † .
If we take p ¼ 0:01, we think there is a 1% probability that a flowset has formed outside the realm of the study.This leads to the relationship where N pre is the number of flowsets formed in the pre-study, N ð0,T is the number of flowsets formed in the study time frame and N is the total number of flowsets.The expected number of flowsets formed in the pre-study is pE M † ½N ¼ 25λ ? .Hence, we can calculate Using Equation ( 4) and substituting the value for λ ?, we find These values of λ and λ ?are fixed for the scoring of the other simulations from the same study.Now we can calculate the loglikelihood for the location and number of lineations component, using A M ðX, TÞ ¼ 35 calculated from our toy example simulation, M, From Figure 2, we can see that flowset formation is only plausible at t ¼ 2 and t ¼ 3. We now calculate the likelihood that the simulation at these two time steps could have formed the flowset, based on the direction of ice flow compared with the lineation direction.We indicate the true location of the observed flowset as x 1 , the corresponding orientation of the flowset as θ 1 and the simulated flow direction at location x 1 and time step t ?as μðx 1 , t ?Þ.At time step 2, we find and at time step 3, Intuitively, this highlights that the ice flow of simulation M at time t ¼ 3 better aligns with the flowset than the ice flow at time t ¼ 2. These values at the two time steps are then summed and multiplied by the rate of formation.The term λ ?2π to account for unrelated flowsets is also added, to calculate the directional component for the likelihood (see Appendix A for more details).
This component can be used in isolation to compare which flowsets contribute the most to the final score and so can indicate which flowsets agree with the simulation the best.
The final log-likelihood for this model simulation is then the difference between the directional and locational components.Note that we need to take the logarithm of the directional component, but the location component is already in the correct form.
Alone, a single value calculated from ( 9) is of limited use.But, through multiple comparisons to simulations, the log-likelihood calculated by LALA can be used to compare ensemble members.Those

| PRACTICAL APPLICATION AND CONSIDERATIONS WHEN USING LALA
Here, we briefly describe a workflow for using LALA.The initial intention and subsequent design of LALA was for comparing multiple ice sheet simulations, to identify fit to the geomorphological record.
However, we anticipate that LALA could also be used to analyse which flowsets fit the best within a single model run, or the timing of best-fit between a simulation and a model.The code for LALA is available from https://github.com/rosiearcher/LALA-Tool.A synthetic example is included with the code.

| Data preparation
Observed data of the flow direction of the former ice sheet in question must be collated.Lineation mapping across the study area should be grouped into flowsets (Clark, 1999): Groups of lineations thought to form at similar times due to their proximity, similar morphology, and orientation (Greenwood & Clark, 2009a;Hughes et al., 2014).

| Code structure
LALA is written in Python (v3.0+) and requires the libraries numpy, pandas, netCDF4, and scipy.The user defined parameters required as inputs are shown in Table 1.As a first step, LALA reads the relevant variables from the input netCDF files (Section 3) and converts them to numpy arrays.Iterating through each time step of the ice sheet simulation, LALA first calculates simulated flow direction.The tool reads in the u and v velocity components from the model simulation, then calculates the angle of the flow using standard trigonometry.To limit calculations to all regions where lineation formation is deemed possible, LALA identifies the locations that exceed the input thickness and velocity thresholds.These are defined by the user as minimum possible conditions for lineation formation (Table 1).The next step is to integrate the area where lineation formation is possible through time for the study area, as well as the area of possible lineation formation for the pre-study.The possible formation area for an ensemble member M † thought to have a plausible size based upon other metrics is calculated and used to adaptively select values for λ and λ ?that will then be applied for the whole ensemble.Next, for each model run, LALA calculates the likelihood of the flowsets forming at each location.For each flowset, times where formation is possible are found and then scored across those time steps to form the directional likelihood (Equation 2).The directional likelihoods are summed over time, for the n flowsets.The final likelihood sums the natural logarithms of each flowset score and sums, to give the log-likelihood, and then the likelihood of locations is subtracted from this to give a score for the simulation (Equation 3).
As output, LALA provides a .xlsxfile with each model simulation number and its associated final log-likelihood score.

| DEFINING THE TOLERANCE BETWEEN MODELLED AND OBSERVED FLOW DIRECTIONS
As demonstrated in the above sections, κ acts as a precision parame- sheet within any grid cell around its overall, more coarsely modelled, mean cell flow direction).We perform an experiment to provide insight into these two components and help us determine suitable values for κ.

| Intermapper comparison experiment
Observations of palaeo-ice flow direction derived from glacial lineations are often mapped manually from digital elevation data (e.g., Finlayson et al., 2014;James et al., 2019;Leger et al., 2020).One form of observational uncertainty when considering palaeo-ice flow direction is that which arises from human error; interpretation of the landform record may vary from person to person.Here, we discuss an experiment performed in collaboration with the geomorphological mapping community to quantify this source of uncertainty (Section 3).
We then use the results from this experiment to examine the uncer- to mappers the selected lineations to map (Figure 5).A summary of the data used, the resolutions, and the number of lineations highlighted are given in Table 2.The participants were asked to draw one crestline for each lineation in the suspected orientation of flow and summary lines for each of the nine boxes that indicated their overall opinion of the ice flow in the specified region.An example of the summary lines of flow direction for one grid cell is shown in Figure 6.We received data from 24 participants.

| Uncertainty from mapping interpretation
We analysed the variability in the reported ice flow direction between the different participants, and for each mapping resolution, across the boxes.We assumed that the summary directions reported by each user within a box followed a von Mises distribution, with each box having a potentially different mean direction μ i but with the same concentration parameter κ 1 shared across all boxes of the same resolution.Using Markov Chain Monte Carlo, we estimated both the mean directions μ i and shared κ 1 to represent the intermapper variability in  3. The concentration (denoting the level of mapping precision) for the lower resolution data was significantly higher than for mapping conducted at a high resolution.This suggests that at a lower resolution, there was a very high level of directional agreement between participants.Several possibilities exist to explain this.This may be a consequence of the nature of the two different study areas (i.e., we happened to choose a much simpler setting for the 30-m test area).Alternatively, the high resolution data  may reveal more detail, representing a more complex landform record that is harder to interpret precisely.Overall however, the values of κ 1 are very high for both resolutions, suggesting that for the areas this experiment was conducted over, the reported summary directions were highly reproducible across mappers and that intermapper variability is not a significant source of uncertainty for model-data comparison.

| Uncertainty between the direction of localised lineations and overall grid cell flow direction
The second source of uncertainty we must consider for our LALA tool is that which results from the comparison of flow directions inferred We selected one mapper's inferred directions for each of the 208 lineations and separated them into 27 groups based on the 5 Â 5 km 2 box that contained them (see Figure 5).We assumed each group/ box could have a different overall mean direction for the lineations, corresponding to the box's overall flow direction.The individual lineations within any box were then modelled by a von Mises distribution, with the appropriate box-dependent mean direction but the same concentration κ 2 for all boxes.MCMC was again performed estimating both the mean direction for each box and the value of κ 2 .The posterior estimate for κ 2 , across all 27 boxes and 208 lineations, was found to be approximately 92.This approach to estimating κ 2 accounts for variability in flow direction within a grid box.For large flowsets, perhaps recording highly variable flow directions, users may wish to either decrease the value of κ 2 to account for such complexity or subset a flowset into multiple grid boxes.

| Combining the two components of directional variability
To use LALA for comparison of observed flowset direction with the output of numerical ice sheet models, we must combine the two components of uncertainty (κ 1 and κ 2 ) described above into a single von flowset direction is very small (with extremely high values of κ 2 ).Consequently, the appropriate overall κ will be dominated by the variability introduced due to the need to compare flow directions at a highly localised lineation scale with the coarser output of an numerical ice sheet model (which only has a resolution on a scale of several km).
We therefore select an overall value of κ ¼ 90 for implemention of LALA.This value will be used us to judge goodness-of-fit between an observed flowset and the flow direction of an ice sheet simulation.In a later version of LALA, we would however endeavour to vary κ spatially, accounting for areas of complex ice flow that likely require a lower concentration value.

| Other potential sources of uncertainty in model-data comparison
The parameter κ captures the uncertainty in the directional component of the likelihood calculated by LALA.However, in any modeldata comparison, there are additional sources of uncertainty, related either to the numerical model utilised or the data itself (Ely et al., 2019).Often, these sources of uncertainty are potentially more difficult to quantify.From a data perspective, an uncertainty remains regarding the mechanics of lineation formation, which are still debated (for a recent summary, see Ely et al., 2022).Though we suggest reasonable conditions under which lineations are formed and preserved, ultimately users of LALA may wish to alter these according to any further insights gained in this field.Furthermore, through the parameter p, LALA accounts for misclassification of lineations from outside of the simulated time period.However, in an ideal situation, such instances of misclassified lineations should be removed from the dataset used to calculate the likelihood of ice sheet simulations.
Model-data comparisons, perhaps facilitated by LALA, may prove a fruitful means for identifying any misfitting data and act as a basis for reinterpreting the lineation record (see Section 3 for an example of this).From a modelling perspective, no ice-flow simulation is a perfect representation of reality.All numerical models abstract the physics of ice flow into approximations, and thus contain structural uncertainty related to the way this is implemented within the model.Thus, the angles of ice flow produced by the model are not perfect representations of the actual ice flow direction.Quantifying the impact of different approximations of ice flow upon the resulting modelled angle is a large task and perhaps a fruitful avenue for future research.A further uncertainty relates to the temporal resolution of the model output.
Again, this is considered by LALA in the parameter p, which considers that lineation forming events could occur between modelled output time steps.However, precisely how modelled angle changes over temporal timescales is unknown, as are the precise timescales over which lineations form.As we gain more insight into these additional sources of uncertainty, future model-data comparison tools or adaptations of LALA may wish to take these into consideration.

| APPLICATION TO THE BRITISH AND IRISH ICE SHEET
To demonstrate the utility of LALA using real model simulations, we here apply LALA to an ensemble of simulations of the British-Irish Ice Sheet (BIIS), from the BRITICE-CHRONO project (Clark et al., 2022).
This 200-member ensemble was run at 5 km horizontal resolution The ensemble experiments varied ten parameters.Our aim here is not to comment on the performance of a particular ice sheet simulation, and its resulting parameter combination, but rather to demonstrate the utility of LALA using a real example of ice flow simulations.
The BRITICE-CHRONO project (Clark et al., 2022) completed ice sheet-wide mapping of glacial landforms of the BIIS (Clark et al., 2018).From the lineation data, flowsets were constructed through grouping lineations thought to indicate similar past flow directions that occurred during the same period of time (Greenwood & Clark, 2009a;Hughes et al., 2014).This flowset data was compiled in GIS, and converted into a netCDF to be used within the LALA tool.
Only flowsets deemed to have formed in the same flow event, referred to as isochronous flowsets, were tested.A total of 94 flowsets were used here, 37 of which were located in England, Scotland and Wales (Hughes et al.2014), and 57 from Northern Ireland and the Republic of Ireland (Greenwood & Clark, 2009a).As the simulated ice sheets we are testing cover the entire late glacial, the time period during which the mapped lineations are thought to have formed, we use the default value of p ¼ 0:01.Thus, in this application of LALA, any observed flowset has a 99% probability that it was formed during the simulated period.
To apply LALA, we must calculate appropriate values for the two parameters which define the rate of lineation formation (Section 3).
For this example, we refer to these as λ ?BC ðxÞ, the rate of lineation formation across the BIIS outside of the simulated period, and λ BC ðx, t, Mðx, tÞÞ, the rate of lineation formation across the BIIS during the study period.For λ ?
BC ðxÞ, we supply LALA with a binary mask defining regions where lineation formation is possible (assigned a value of 1, which LALA later translates to λ ? BC [Section 3]), and impossible (assigned a value of 0).The latter ruled-out regions were defined as those where we deemed there was insufficient sediment for lineation formation, offshore regions as these were not covered by Greenwood and Clark (2009a) and Hughes et al. (2014), and regions beyond 50 km of the presumed limit of the BIIS (Clark et al.2022).This provided us with a total number of grid cells where lineations might occur AðX Þ ¼ 7801.Using this, and the number of observed flowsets (n ¼ 94), we can then calculate λ ?BC as To estimate the parameter λ BC , we need to determine AðX, TÞ, the time-integrated area which meets the lineation formation criteria we previously formulated.One could envisage estimating AðX, TÞ via a number of means, dependent upon the user's opinion regarding lineation formation and the availability of information regarding the overall extent and timing of glaciation across an area.For example, an empirical reconstruction could be used to inform about the time for which each cell was covered by ice.However, such reconstructions often lack information about the prevalence of other conditions thought to be important for lineation formation, such as basal thermal regime.
Our approach here is to utilise an initial simulation, M † , that scored well at matching the independent data of ice margin position (Li et al., 2008) and geochronological constraints (Ely et al., 2019).M † is therefore of the approximate size and produced reasonable ice coverage for the BIIS.For simulation M † , the total number of grid cells (area) These two rates, λ ?BC ðxÞ and λ BC ðx, t,Mðx, tÞÞ, were then used in LALA to calculate the log-likelihood of each simulation from the ensemble.For the remaining input parameter, κ, we use the value 90, derived from our observational uncertainty experiment (Section 3).
To assess the utility of LALA, we visually inspected how a range of simulations compared with the reconstructed ice extents of Clark et al. (2022), at the 1-ka temporal resolution at which the reconstructed limits are provided.An example of this comparison, for the simulations that had the highest and lowest log-likelihoods, is shown in Figure 7.Note that for illustrative purposes in Figure 7, we demonstrate this for a single time step, 21-ka BP, whereas our visual assessment accounted for the time variance of the ice sheet geometry.Visual assessment of all 200 simulations within the ensemble showed that the best-fitting simulation identified by LALA performed reasonably well in comparison to Clark et al. (2022).Furthermore, the worst simulation exceeds the maximum reconstruction almost all the way around the extent, and is comparable in performance to other poorly performing models identified through qualitative means.The best performing simulation (Figure 7a) has an extent that matches the reconstructed ice margin extent for the majority of its perimeter.In comparison, the extent of the worst performing model (Figure 7b) is consistently too large when compared to the empirical reconstruction.This conformity between reconstructed and simulated ice sheet geometry reflects how LALA produces a higher log-likelihood for simulations that produce lineation forming conditions (velocity, ice thickness, and grounded extent) over the model domain, and penalises simulations that spread over an extent which is too large.By way of intuitively noting how well each simulation matches purely the flow direction of flowsets, we also isolated only the component of the loglikelihood related to direction and summed this on a per-simulation basis (Section 3).The best fitting simulation (Figure 7a) also showed the overall best directional score, and vice versa for the worst fitting simulation.As LALA considers more than just flow direction alone, a quantitative comparison between LALA and AFDA is not currently possible.However, we note that higher-ranked simulations by LALA were also those assessed to fit more flowsets by AFDA.
A further utility of LALA is to consider which flowsets are more readily matched, and those which are less often matched.To do so, we isolate only the component of the log-likelihood that relating to flow direction (Equation 2).When scoring a whole model simulation, LALA sums the values obtained by Equation ( 2) to calculate the final log-likelihood.To give an indication of the best matched flowsets across the ensemble, we instead keep the Equation ( 2) values separate.Then, each flowset will have a value for each model simulation.
We summed these values from each of the simulation on a perflowset basis.The nature of this comparison may provoke iterative conversations between those who collect and collate data, and the numerical modeller: is the match, or lack of, a consequence of the model, the nature of the of data, or some other factor?In Figure 8, we show the most frequently matched flowsets, and those that were matched least frequently.Although we do not revisit the evidence for these flowsets here, there are some notable features of this map.

| SUMMARY
Rigorous tools for comparing numerical simulations of palaeo-ice sheet behaviour and available data constraints can greatly aid in model improvement.Here, we present a new method for determining the log-likelihood of an ice sheet simulation given a set of mapped subglacial lineations.This tool, which we name the Likelihood of Accordant Lineations Analysis (LALA) tool, considers both the orientation and spatial distribution of lineations across a palaeo-ice sheet area.To achieve this, we used the statistical underpinning of considering lineation formation to be a marked, inhomogeneous, Poisson point process.Unlike previous tools, LALA provides a continuous numerical output, suited for emulating ensembles of simulations.The main purpose of LALA is to provide a means for comparing ice sheet simulations within an ensemble, assigning each simulation with a log-likelihood, such that the best-fitting simulations can be identified.To demonstrate the underpinning of LALA, we provide a toy example of its application.
Several parameters are required as a user input to run LALA.As LALA is designed to be adaptable, these can be changed by the user depending upon assumptions of lineation formation and to adapt to different experimental circumstances.Here, we provide a number of reasonable starting points for assigning these parameters.This includes an analysis of an experiment which quantified the observational error that arises from multiple expert mappers interpreting former ice flow directions from the landform record.We find this uncertainty to be small, especially compared to the uncertainty that arises from the differing resolutions of subglacial lineations (100s of m) to the scale of an ice sheet model grid (several km), which we find to be much larger.Finally, to demonstrate the utility of LALA, we apply the tool to an ensemble of simulations of the British-Irish Ice Sheet and observations of lineations from the literature.This application highlights how LALA can be used to distinguish between simulations within an ensemble, and in the future find more plausible parameter spaces for ice sheet simulations.Conditional on lineation formation occurring at time t and location x, a flowset will be created that has a direction related to the basal flow direction of the ice sheet at point of creation.The flowsets therefore constitute a marked Poisson point process, with the mark being their direction.Specifically, we assume that ϕ, the direction of the flowset, will be distributed (independently of other flowsets) according to a von Mises distribution: Here κ 1 denotes the concentration of the true flowset around the ice flow direction.
Unrelated flowsets poisson point process (to account for lineations forming outside of the study region and time period) To allow for the possibility that some of the observed flowsets forming outside of the study we introduce an additional, independent, Poisson point process over the same region X but over a (pseudo)time period ðÀ1,0 with intensity λ ?ðxÞ.Here the intensity depends only upon the location x, typically, we would expect λ ?ðxÞ ¼ 0 in ocean locations, or where the topography is such that flowsets are impossible.We assume that formation of lineations created from this additional point process lead to flowsets that lie uniformly in any direction, that is, the direction ϕ of any flowset created in this pseudo-time period has probability density This additional Poisson point process aims to represent the possibility a flowset may have been created outside the realm of the study (before the time t ¼ 0 when our modelling starts).In any subregion A, the number of (unrelated) flowsets also follows a Poisson distribution with mean

Observational model
We assume that in our region under study X , we observe n flowsets in locations x 1 , …, x 1 .Note that the region X includes only the areas mapped for lineations (regions where the presence/absence of lineations is unknown do not form part of X).Each flowset has an estimated direction θ i based on mapping.We assume that the observed θ i can either relate to the underlying angle ϕ i of the flowset, or the opposite direction ϕ i þ π.The scientist describing the flowset may equally report the direction as being upstream or downstream due to their morphology.We also assume that there is also some measurement uncertainty in the estimate θ i so that the density of θ i jϕ i is Here, κ 2 relates to the precision of the observed mapping (i.e., how precise a user is at reporting the true lineation direction).
For simplicity in modelling, we combine the observational uncertainty in lineation direction reported by the mapper, with the variability around the ice flow direction at time of formation.When the lineation is formed by an ice sheet M, a slight approximation leads to a final observation model for θ of f M ðθ i jμðx i , t ?i ÞÞ ¼ 0:5f vonM ðθ i jμðx i , t ?i Þ, κÞ þ 0:5f vonM ðθ i jμðx i , t ?i Þ þ π, κÞ, ð †Þ where t ?i denotes the time of the formation that led to the flowset at x i and κ is a concentration term combining both the observational and regridding components of variability.When the lineation is formed outside of the study, we still have f out ðθÞ ¼ 1 2π : Log-likelihood of fx i , θ i g n i¼1 Due to their independence, we can combine the point process for ice sheet M (over times T ¼ ð0,T) and the outlier point process (to represent potential earlier ice sheets) to obtain another marked inhomogenous Poisson point process.For this combined process, the likelihood of seeing the n flowsets in locations g e ÀðΛðX ,ð0,TÞþΛ ?ðX ÞÞ , where Λðx i Þ ¼ Ð T 0 λðx i , t, Mðx i , tÞÞdt.Conditional on the flowset locations, we can also calculate the likelihood of the observed directions θ i for i ¼ 1, …, n.This is achieved by conditioning on t ?
i the time of the flowset formation at location x i .We find fort ?ð0, T, i:e:; formed by an ice sheetM; fort ?ðÀ1,0, i:e:; unrelated to the study: Conditional on the time t ?i that the flowset in location x i is formed, the distribution of θ i is known.Hence, we find the likelihood of θ i jx i ,M to be where f M ðθ i jμðx i , t ?i ÞÞ is given in (A1).Putting together both the location and direction information, we find that the final log-likelihood for ice sheet M is

Discretisation
In practice, the ice sheet model simulation will consist of individual time steps defined on a grid rather than over continuous time and space.Hence when testing an ice sheet model simulation, the integrals above (νðθ i jx i , MÞ, ΛðX, ð0, TÞ and Λ ?ðXÞ) will be replaced by discrete sums (over the grid cells, time steps, or both).

Selecting probability an observed lineation is unrelated
We suggest that, rather than choose a value for λ ?directly, we instead specify p, our prior probability that an observed flowset arises from outside the study's scope.The number of flowsets created by the ice sheet Mðx, tÞ over the study period T and region X follows a Poisson distribution with mean ΛðX , ð0, TÞ.The number of flowsets created by our (pseudo-time) ice sheet to represent formation for reasons not considered in the study has a Poisson distribution with mean Λ ?ðXÞ.We can therefore obtain our chose value of p by setting ΛðX, ð0, TÞ: Typically, we choose λ ?ðx i Þ ¼ λ ?for regions where lineations are possible, and zero otherwise.Hence, Λ ?ðXÞ ¼ λ ?AðXÞ where AðXÞ is the area where lineations could form within the overall study region X.Our default is to select p ¼ 0:01.
Choosing λðx i , t ?i , Mðx i , t ?i ÞÞ In our use of LALA, we have selected λðx i , t ?i , Mðx i , t ?i ÞÞ ¼ 0 whenever the location x i is not covered by grounded ice that is at least 10 m thick and travelling at 10 ms À1 or higher.Effectively, it is not possible for the ice sheet to form lineations. Otherwise, when the ice sheet has the potential to form lineations at location x, we have chosen it to be spatially and temporally uniform, that is, λðx i , t ?i , Mðx i ,t ?i ÞÞ ¼ λ.If we select such a uniform λ, then ΛðX, ð0, TÞ ¼ λA M ðX , TÞ where 2πI 0 ðκÞ is a von Mises distribution.Here, κ is a measure of concentration (i.e., how close the directions of the lineations lie to the underlying ice sheet flow) and we have measured angles in radians (where 360 is 2π radians).Examples of the von Mises probability density function for different values of κ are shown in Figure 1.Large values of κ mean that the lineations created during formation will lie tightly around the direction of ice flow, while smaller values of κ permit greater variability in the direction of the lineations.We discuss the selection of κ in Section 3.
The probability density function of three various von Mises distributions with the same mean, 0, and three different shape parameter values, κ.As κ increases, then the distribution becomes more concentrated around its mean.Here, angles are shown in degrees (360 is 2π radians).[Color figure can be viewed at wileyonlinelibrary.com] formation are at the discretion of the user; however, the stated conditions are used in this study.In such a two-state model, we are however still required to select an appropriate value for λ during times of potential lineation formation.We propose a sensible selection obtained by calculating, for an initial ice sheet simulation M † that approximately accords with other palaeo-information such as ice extent or volume, the expected number of flowsets E M † ½N according to our Poisson process model.This

κ
is the concentration of the observed lineations around the ice sheet flow direction during lineation formation.As shown in Figure1, the value of κ changes the shape of the von Mises distribution, demanding a higher level of directional accordance if κ is large.The appropriate value of κ, a measure of the acceptable direction for an ice sheet simulation, is dependent upon two components.The first component relates to the variation in observer measurement of the true lineation direction-which any mapper introduces a level of individual measurement error.The second pertains to the uncertainty regarding the relationship between the true (mean) direction of a flowset and that of the numerically modelled ice flow of simulation M within the grid cell.This second component itself consists of two parts-which we are gridding (summarising) multiple lineations orientations into a single flowset at the scale of an ice sheet grid and that the numerical ice sheet model output is somewhat coarse and may not represent localised variation in ice flow direction within a grid cell (within any grid cell, the direction of ice flow at a more localised scale may vary from the overall flow direction within the cell).We consider these two components and practical means for choosing a value of κ in Section 3.
, we show the grid cells at each of the three time steps that could feasibly form a flowset, following the conditions we have previously set out regarding grounded ice, minimum velocity, and thickness.The location of the flowset is indicated in Figure 2 by a cross.Figure 3 indicates the difference in simulated flow direction and the observed direction inferred from the geomorphological evidence.The observed direction for flow, assumed to be parallel to past flow direction is taken to be 45 .While the directions are shown in degrees for this illustrative example, calculations in LALA (and all of the equations presented in this paper) are formed a flowset, and so AðX Þ ¼ 25.We use a value of A M † ðX, TÞ ¼ 42 for the time-integrated plausible area for lineation formation for the simulation M † .The value for A M † ðX , TÞ has been chosen arbitrarily for F I G U R E 2 Three simulated time steps created to illustrate the use of plausible areas for lineation formation.The grid represents a 5 Â 5 domain and each box is coloured either blue or pink dependent on whether the location is deemed to be plausible or implausible, respectively, for the formation of lineations.The location of the flowset to be scored is marked with a black cross.[Color figure can be viewed at wileyonlinelibrary.com] the toy example.A more formal method for choosing A M † ðX ,TÞ is discussed in Section 3. The values of λ and λ ?are the same for each model simulation in the ensemble, using the fixed model simulation M † .Under our Poisson process model, the number of flowsets we would expect to see with simulation M † over our three time step study period, combined with the prestudy, is 1 as there is one flowset under consideration.In total therefore, with simulation M † we would

F
I G U R E 3 Three time steps showing the difference between the simulated ice flow direction and the inferred lineation direction created to illustrate how LALA scores an ice sheet simulation.The observed flowset is represented by the blue boxes with the black arrows showing the orientation given at the middle of the flowset.The simulated flow direction is indicated by the grey lines.The difference between the observed and modelled flow direction is indicated for each time step.[Color figure can be viewed at wileyonlinelibrary.com] simulations with greater log-likelihoods would be considered more likely to have generated the observed lineations, while those with lower log-likelihoods are less likely.
Currently, a single pixel located at the centre of each flowset should be used.Expert interpretation is required for choosing this point for complicated flowsets and special attention should be given to the relative size of the flowsets compared with the model resolution.The flowset data should take the form of a netCDF grid, with the same resolution and extent as the model simulation output.Each flowset should be represented by one layer of the netCDF file, making the data a three-dimensional array with the flowset number as a dimension.This makes the individual flowsets more accessible and removes any overlapping lineation data (i.e., where lineations from two flowsets occupy the same grid cell, perhaps cross-cutting).A netCDF of the ice sheet simulation output is also required.This should contain the following variables, with their standard names in PISM provided in brackets: basal velocity magnitudes in the u and v directions (uvel and vvel), ice thickness (thk), ice speed (velsurf_mag), and a mask of ice extent (mask).Finally, users can specify a third input netCDF file, which contains information pertaining to the regions where lineations are likely to have formed.In the absence of this file, the whole domain is presumed equally likely to have conditions conducive to lineation formation in the prestudy.
ter for comparing simulated and observed flow directions.Higher values of κ will produce a lower tolerance for model-data misalignment, and vice versa.There are two aspects of the model-data comparison procedure that contribute to the value of κ that the user should define when using LALA: (i) Observational variability in the reported mean orientation of a flowset; and (ii) uncertainty in the relationship between the mean flowset orientation and the (typically coarse) output of an ice sheet simulation (consisting of both the uncertainty introduced when regridding lineation data to model resolution, and the additional localised variability in the flow of an ice T A B L E 1 Inputs that need to be provided by the user to score ice sheet model simulations using LALA.
Figure 4.The areas were each split into nine equally sized 5 Â 5 km 2 boxes.A selection of lineations were starred in each area to indicate Overview of areas used in the lineation mapping experiment.[Color figure can be viewed at wileyonlinelibrary.com] assessing flowset direction for a particular resolution.Further details on this method can be found in the Appendix.We refer to this experiment as the Multiple Means experiment.For the Multiple Means experiment, we first considered the two high resolution (2 m) mapping results separately, later treating them as a single population, which allowed us to get a single value of κ 1 representing this resolution.Results are shown in Table

F
I G U R E 5 Maps of the three areas given to the observational error experiment participants.The black points show the location of the lineations participants were asked to map.The 3 Â 3 grid shows the sections where summary lines were requested.(a) Area 1 located in northern Finland.(b) Area 2 located in central Finland.(c) Area 3 in western Russia.[Color figure can be viewed at wileyonlinelibrary.com]T A B L E 2 Data used in the observational error experiment.
at a lineation scale (generally 100s of m in scale) with those directions provided by an ice sheet simulation that typically has a much coarser output, representing mean ice flow direction at a several km scale.In the case of the model simulations from the BRITICE-CHRONO project (Clark et al., 2022).The simulations are transient and run from 31 to 15 ka, outputting every 100 years, with model domain covering Britain and Ireland.To assess the uncertainty in model-data comparison this reduction introduces, we considered the variability of the individual lineations within a 5 Â 5 km 2 box around their overall mean direction.Boxes of this size were chosen to match the resolution of the BRITICE-CHRONO ensemble.
Mises concentation parameter κ (see Appendix B).Our experiment above indicates that intermapper (observer) variability in the reported F I G U R E 6 One of the nine boxes from area 2 where the mapping participants seemed to identify two overlapping directions.Left: Hillshade of the area requested to be mapped, with the specific lineations to be mapped starred.Right: The summary lines from the 24 participants are indicated in blue.The extreme angles are highlighted in purple to emphasise the difference between mappers.As the lines here indicate the inferred summary orientation for each box, and not the specific marked lineations, the lines are placed centrally within the box.[Color figure can be viewed at wileyonlinelibrary.com]T A B L E 3 Values obtained from MCMC analysis of reported summary directions of lineations across boxes within our intermapper experiment.Results are presented separately for three different regions, two at a 2-m resolution and one at a 30-m resolution; as well as combining the two 2-m resolution areas.Each box was permitted to have a different overall mean direction.
using the Parallel Ice Sheet Model (PISM) (Winkelmann et al., 2011)a hybrid shallow ice and shallow shelf approximation numerical ice sheet model.The simulations are transient and run from 31 to 15 ka with model domain covering Britain and Ireland(Clark et al., 2022).

First
, a number of the poorly fitted flowsets are narrow and small features (e.g., FS58 in the Republic of Ireland, FS24 on the border between Northern Ireland and the Republic of Ireland, and FS47 in Northern Ireland).These flowsets were perhaps poorly chosen for model-data comparison, as their narrow size may be below what is feasibly resolved by the model.In western Scotland, FS24 is composed of two distinct patches of lineations, which cross the fjord topography.Perhaps these lineations have been incorrectly grouped together, or the model is unable to recreate flow counter to the valley topography.The reasons for the well matched flowsets may be highly variable, as their character and distribution seem to lack geographic, topographic, or geometric similarity (Figure 8).Perhaps these flowsets F I G U R E 7 Modelled extents of the BIIS at 21 ka BP in two different simulations.The final log-likelihood score for each simulation is indicated below.Note that this figure shows one time step, but the final scores are cumulative over all time steps.The lines denote the maximum (red), optimum (yellow) and minimum (green) reconstructed extent of the BIIS from Clark et al. (2022).The faint white-line is the grounding line or ice margin.(a) The best-fit simulation.Note how much of the perimeter of the ice sheet is close to the 'optimum' reconstruction.(b) The worstfitting simulation.Note how the ice sheet extent is too large compared to the empirical reconstruction.[Color figure can be viewed at wileyonlinelibrary.com] are looser constraints on ice sheet models, or the ensemble simulation is particularly good at replicating these flow directions.Such possibilities should be considered by future work.

F
I G U R E 8 The highest scoring flowsets across all model simulations are indicated in blue.The lowest scoring flowsets shown in yellow.The timeintegrated score for each simulation given by LALA is below.[Color figure can be viewed at wileyonlinelibrary.com] https://doi.org/10.1002/esp.5658APP E NDIX A: LOG-LIKELIHOOD OF A POISSON POINT PROCESS Note that all angles θ described below are measured in radians rather than degrees; 360 is equal to 2π radians.An ice sheet point process Assume that, for a given ice sheet Mðx, tÞ over a study period T ¼ ð0, T and region X, lineation formation follows an inhomogeneous Poisson point process with intensity function λ M ¼ λðx, t,Mðx, tÞÞ.For any subregion A and time interval ða, b.The number flowsets formed NðA, ða,bÞ will follow a Poisson distribution with mean ΛðA, ða, bÞ ¼

A
M ðX , TÞ effectively integrates (over time and space) the region when the ice sheet M can form lineations, AðX , TÞ ¼ ðð Lineation formation possible with ice sheet M dtdx:The total expected number of flowsets (once we include the pseudo-time sheet) is Nðλ, λ ?, MÞ ¼ λA M ðX , TÞ þ λ ?AðXÞ.To choose a suitable λ we first select a value of p.This will reduce Nðλ, λ ?, MÞ to a dependence upon only λ.We calculate this expected value Nðλ, MÞ for a single selected ice sheet M ?(ideally one which has been identified as plausible based on other palaeo-data) and match this to n the observed number of flowsets.AP PE NDIXB: USING Markov chain Monte Carlo (MCMC) TO ESTIMATE THE APPROPRIATE VON MISES κ PARAMETERS TO DESCRIBE VARIATION IN FLOWSET DIRECTION We used the MCMC to estimate a set of von Mises parameters that adequately represents the flowset variation observed with our experimental mapping data in Section 3.These parameters are then used as the default choice for LALA.A Markov Chain, a path which moves through different states, was generated according to a Metropolis-Hastings (MH) algorithm.The MH algorithm sampled a random point from a proposed density function.The point is then either accepted as the next state in the chain, or rejected, and a new sample point is tested.This is repeated until the distribution of acceptable points stabilises to form a posterior distribution.From this sampled posterior distribution, metrics such as the mean and variance can be calculated.Several thousands of iterations are performed until the Markov Chain reaches its stationary distribution.Examples of the Markov chains used for our analysis of area 2 are shown below in Figure B1.F I G U R B Sample chains from the MCMC for the 2 m resolution data in area 2. Left: Considers there to be one mean across all the boxes.Right: Considers a separate mean for each box, but the same precision.The red line shows the mean that the posterior distribution tends towards.[Color figure can be viewed at wileyonlinelibrary.com] The values of λðx, t, Mðx, tÞÞ and λ ?ðxÞ can also be changed.For λðx, t,Mðx, tÞÞ especially, the conditions under which the user believes lineations are able to form can be altered.If the assumptions include the same properties as we look at here, that is, thickness, velocity, and grounded ice, but different values are wanted, there is a predefined function that can easily change the values.If extra conditions are to be included, the easiest method is to create a binary map with the area at each time step that meets the user's extended conditions and use this in place of the area calculated in the LALA tool.
† ðX, TÞ ¼ 6094155.We can use this to select λ BC for our ensemble