diff --git a/Docs/ELA-model_tex/Figures/ELAs.png b/Docs/ELA-model_tex/Figures/ELAs.png
new file mode 100644
index 0000000..b58df24
Binary files /dev/null and b/Docs/ELA-model_tex/Figures/ELAs.png differ
diff --git a/Docs/ELA-model_tex/Figures/Hypsometry.png b/Docs/ELA-model_tex/Figures/Hypsometry.png
new file mode 100644
index 0000000..ca05d2a
Binary files /dev/null and b/Docs/ELA-model_tex/Figures/Hypsometry.png differ
diff --git a/Docs/ELA-model_tex/Figures/Ice.png b/Docs/ELA-model_tex/Figures/Ice.png
new file mode 100644
index 0000000..5bd9917
Binary files /dev/null and b/Docs/ELA-model_tex/Figures/Ice.png differ
diff --git a/Docs/ELA-model_tex/Figures/Width.png b/Docs/ELA-model_tex/Figures/Width.png
new file mode 100644
index 0000000..8467407
Binary files /dev/null and b/Docs/ELA-model_tex/Figures/Width.png differ
diff --git a/Docs/ELA-model_tex/Figures/arc-flow.png b/Docs/ELA-model_tex/Figures/arc-flow.png
new file mode 100644
index 0000000..646ca17
Binary files /dev/null and b/Docs/ELA-model_tex/Figures/arc-flow.png differ
diff --git a/Docs/ELA-model_tex/Figures/graph-abstract.png b/Docs/ELA-model_tex/Figures/graph-abstract.png
new file mode 100644
index 0000000..40f731b
Binary files /dev/null and b/Docs/ELA-model_tex/Figures/graph-abstract.png differ
diff --git a/Docs/ELA-model_tex/Figures/mermaid-figs.md b/Docs/ELA-model_tex/Figures/mermaid-figs.md
new file mode 100644
index 0000000..1a3cbdc
--- /dev/null
+++ b/Docs/ELA-model_tex/Figures/mermaid-figs.md
@@ -0,0 +1,101 @@
+# Mermaid diagrams
+
+A space to keep archival code to generate mermaid diagrams for use in the paper.
+Go to the [live editor](https://mermaid-js.github.io/mermaid-live-editor/) in order to generate .png of the charts.
+
+## Graphical abstract
+
+```mermaid
+graph LR
+ subgraph Inputs
+ a1>Bed topography
measurements]
+ a2>Glacier width
measurements]
+ a3>Model parameters]
+ end
+ subgraph Intermediate Components
+ b1{{Modeled bed
topography}}
+ b2{{Modeled ice
thickness}}
+ b3{{Modeled width}}
+ end
+ subgraph Model
+ c{ELA_calc.m}
+ end
+ subgraph Outputs
+ d1([ELA estimate])
+ d2([ELA error])
+ end
+ a1 ==>|+Error| b1
+ a3 ==> b1
+ b1 ==>|+Error| b2
+ a3 ==> b2
+ b2 ==> b3
+ a2 ==>|+Error| b3
+ a3 ==> b3
+ b1 ==>|+Error| c
+ b2 ==>|+Error| c
+ b3 ==>|+Error| c
+ c ==> d1
+ c ==>|+Monte Carlo| d2
+
+classDef default fill:#5cfcfc,stroke:#333,stroke-width:1px
+```
+
+```json
+{
+ "theme": "forest",
+ "flowchart": {
+ "htmlLabels": false,
+ "useMaxWidth": false,
+ "rankSpacing": 15,
+ "nodeSpacing": 100
+ }
+}
+```
+
+## ArcGIS workflow
+
+```mermaid
+graph LR
+ subgraph Inputs
+ iA[Aerial
imagery]
+ iB[DEM]
+ end
+ A["Glacier
centerline
(freehand)"]
+ B[Glacier
outline]
+ iA --> A
+ iA --> B
+ A --> tA
+ subgraph Tempory .xls file
+ tA[Export
transect
elevation]
+ tB[Export
transect
ice surface]
+ end
+ iB --> tA
+ D["Find ice surface
elevation"]
+ tA --> D
+ D --> tB
+ C[Orthogonal
centerline
transect]
+ tB --> C
+ B --> E[Clip width to
outline/
topography]
+ C --> E
+ iB --> E
+ subgraph Output
+ oA[Export to .csv]
+ end
+ E --> oA
+
+classDef default fill:orange,stroke-width:2.5px
+classDef matlab fill:#5cb1fc, stroke-width:2.5px
+class D matlab
+```
+
+```json
+{
+ "theme": "forest",
+ "flowchart": {
+ "htmlLabels": false,
+ "useMaxWidth": false,
+ "rankSpacing": 15,
+ "nodeSpacing": 50
+ }
+}
+```
diff --git a/Docs/ELA-model_tex/citations.bib b/Docs/ELA-model_tex/citations.bib
new file mode 100644
index 0000000..32980fb
--- /dev/null
+++ b/Docs/ELA-model_tex/citations.bib
@@ -0,0 +1,287 @@
+
+@article{benn_mass_2000,
+ title = {Mass balance and equilibrium-line altitudes of glaciers in high-mountain environments},
+ volume = {65-66},
+ issn = {1040-6182},
+ url = {http://www.sciencedirect.com/science/article/pii/S1040618299000348},
+ doi = {10.1016/S1040-6182(99)00034-8},
+ abstract = {The mass-balance characteristics of glaciers in high-mountain environments complicate the relationship between glacier equilibrium-line altitudes (ELAs) and climatic variables such as precipitation and air temperature. Therefore, methods of ELA reconstruction employed in low-relief environments are commonly not applicable in high mountains, or require some modification. We review the concept of the ELA, with reference to the mass balance of a range of glacier types found in high-mountain regions. We examine the applicability of several commonly used methods of ELA reconstruction for different glacier types, and propose some general principles to guide the choice of appropriate methods.},
+ language = {en},
+ urldate = {2020-04-28},
+ journal = {Quaternary International},
+ author = {Benn, Douglas I and Lehmkuhl, Frank},
+ month = apr,
+ year = {2000},
+ pages = {15--29}
+}
+
+@article{farinotti_method_2009,
+ title = {A method to estimate the ice volume and ice-thickness distribution of alpine glaciers},
+ volume = {55},
+ issn = {0022-1430, 1727-5652},
+ url = {https://www.cambridge.org/core/journals/journal-of-glaciology/article/method-to-estimate-the-ice-volume-and-icethickness-distribution-of-alpine-glaciers/40D2B92744F6512E8DEC142DE6BEA471},
+ doi = {10.3189/002214309788816759},
+ abstract = {Sound knowledge of the ice volume and ice-thickness distribution of a glacier is essential for many glaciological applications. However, direct measurements of ice thickness are laborious, not feasible everywhere and necessarily restricted to a small number of glaciers. In this paper, we present a method to estimate the ice-thickness distribution and the total ice volume of alpine glaciers. This method is based on glacier mass turnover and principles of ice-flow mechanics. The required input data are the glacier surface topography, the glacier outline and a set of borders delineating different ‘ice-flow catchments’. Three parameters describe the distribution of the ‘apparent mass balance’, which is defined as the difference between the glacier surface mass balance and the rate of ice-thickness change, and two parameters define the ice-flow dynamics. The method was developed and validated on four alpine glaciers located in Switzerland, for which the bedrock topography is partially known from radio-echo soundings. The ice thickness along 82 cross-profiles can be reproduced with an average deviation of about 25\% between the calculated and the measured ice thickness. The cross-sectional areas differ by less than 20\% on average. This shows the potential of the method for estimating the ice-thickness distribution of alpine glaciers without the use of direct measurements.},
+ language = {en},
+ number = {191},
+ urldate = {2020-04-28},
+ journal = {Journal of Glaciology},
+ author = {Farinotti, Daniel and Huss, Matthias and Bauder, Andreas and Funk, Martin and Truffer, Martin},
+ year = {2009},
+ note = {Publisher: Cambridge University Press},
+ pages = {422--430}
+}
+
+@article{haeberli_application_1995,
+ title = {Application of inventory data for estimating characteristics of and regional climate-change effects on mountain glaciers: a pilot study with the {European} {Alps}},
+ volume = {21},
+ issn = {0260-3055, 1727-5644},
+ shorttitle = {Application of inventory data for estimating characteristics of and regional climate-change effects on mountain glaciers},
+ url = {https://www.cambridge.org/core/journals/annals-of-glaciology/article/application-of-inventory-data-for-estimating-characteristics-of-and-regional-climatechange-effects-on-mountain-glaciers-a-pilot-study-with-the-european-alps/2F94393AD2866405FF3490C74D4EAE2D},
+ doi = {10.3189/S0260305500015834},
+ abstract = {A parameterization scheme using simple algorithms for unmeasured glaciers is being applied to glacier inventory data to estimate the basic glaciological characteristics of the inventoried ice bodies and simulate potential climate-change effects on mountain glaciers. For past and potential climate scenarios, glacier changes for assumed mass-balance changes are calculated as step functions between steady-state conditions for time intervals that approximately correspond to the characteristic dynamic response time (a few decades) of the glaciers. In order to test the procedure, a pilot study was carried out in the European Alps where detailed glacier inventories had been compiled around the mid-1970s. Total glacier volume in the Alps is estimated at about 130 km3 for the mid-1970s; strongly negative mass balances are likely to have caused a loss of about 10–20\% of this total volume during the decade 1980–90. Backward calculation of glacier-length changes using a mean annual mass balance of 0.25m w.e.a−1 since the end of the “Little Ice Age” around 1850 AD gives considerable scatter but satisfactory overall results as compared with long-term observations. The total loss of Alpine surface ice mass since 1850 can be estimated at about half the original value. An acceleration of this development, with annual mass losses of around 1 m a−1 or more as anticipated from IPCC scenario A for the coming century, could eliminate major parts of the presently existing Alpine ice volume within decades.},
+ language = {en},
+ urldate = {2020-04-28},
+ journal = {Annals of Glaciology},
+ author = {Haeberli, Wilfried and Hoelzle, Martin},
+ year = {1995},
+ note = {Publisher: Cambridge University Press},
+ pages = {206--212},
+ file = {Haeberli_Hoelzle_1995_Application of inventory data for estimating characteristics of and regional.pdf:/home/durbank/ODrive/Zotero attachments/Haeberli_Hoelzle_1995_Application of inventory data for estimating characteristics of and regional.pdf:application/pdf}
+}
+
+@article{more_computing_1983,
+ title = {Computing a {Trust} {Region} {Step}},
+ volume = {4},
+ issn = {0196-5204},
+ url = {https://doi.org/10.1137/0904038},
+ doi = {10.1137/0904038},
+ abstract = {We propose an algorithm for the problem of minimizing a quadratic function subject to an ellipsoidal constraint and show that this algorithm is guaranteed to produce a nearly optimal solution in a finite number of iterations. We also consider the use of this algorithm in a trust region Newton's method. In particular, we prove that under reasonable assumptions the sequence generated by Newton's method has a limit point which satisfies the first and second order necessary conditions for a minimizer of the objective function. Numerical results for GQTPAR, which is a Fortran implementaton of our algorithm, show that GQTPAR is quite successful in a trust region method. In our tests a call to GQTPAR only required 1.6 iterations on the average.},
+ number = {3},
+ urldate = {2020-04-28},
+ journal = {SIAM Journal on Scientific and Statistical Computing},
+ author = {Moré, Jorge J. and Sorensen, D. C.},
+ month = sep,
+ year = {1983},
+ keywords = {ellipsoidal constraint, global convergence, Newton's method, trust region},
+ pages = {553--572},
+ file = {Moré_Sorensen_1983_Computing a Trust Region Step.pdf:/home/durbank/ODrive/Zotero attachments/Moré_Sorensen_1983_Computing a Trust Region Step.pdf:application/pdf}
+}
+
+@article{machguth_exploring_2008,
+ title = {Exploring uncertainty in glacier mass balance modelling with {Monte} {Carlo} simulation},
+ volume = {2},
+ issn = {1994-0416},
+ url = {https://www.the-cryosphere.net/2/191/2008/},
+ doi = {https://doi.org/10.5194/tc-2-191-2008},
+ abstract = {{\textless}p{\textgreater}{\textless}strong{\textgreater}Abstract.{\textless}/strong{\textgreater} By means of Monte Carlo simulations we calculated uncertainty in modelled cumulative mass balance over 400 days at one particular point on the tongue of Morteratsch Glacier, Switzerland, using a glacier energy balance model of intermediate complexity. Before uncertainty assessment, the model was tuned to observed mass balance for the investigated time period and its robustness was tested by comparing observed and modelled mass balance over 11 years, yielding very small deviations. Both systematic and random uncertainties are assigned to twelve input parameters and their respective values estimated from the literature or from available meteorological data sets. The calculated overall uncertainty in the model output is dominated by systematic errors and amounts to 0.7 m w.e. or approximately 10\% of total melt over the investigated time span. In order to provide a first order estimate on variability in uncertainty depending on the quality of input data, we conducted a further experiment, calculating overall uncertainty for different levels of uncertainty in measured global radiation and air temperature. Our results show that the output of a well calibrated model is subject to considerable uncertainties, in particular when applied for extrapolation in time and space where systematic errors are likely to be an important issue.{\textless}/p{\textgreater}},
+ language = {English},
+ number = {2},
+ urldate = {2020-04-28},
+ journal = {The Cryosphere},
+ author = {Machguth, H. and Purves, R. S. and Oerlemans, J. and Hoelzle, M. and Paul, F.},
+ month = dec,
+ year = {2008},
+ note = {Publisher: Copernicus GmbH},
+ pages = {191--204},
+ file = {Machguth et al_2008_Exploring uncertainty in glacier mass balance modelling with Monte Carlo.pdf:/home/durbank/ODrive/Zotero attachments/Machguth et al_2008_Exploring uncertainty in glacier mass balance modelling with Monte Carlo.pdf:application/pdf}
+}
+
+@book{oerlemans_minimal_2008,
+ address = {Utrecht},
+ title = {Minimal glacier models},
+ isbn = {978-90-6701-022-1},
+ language = {English},
+ publisher = {Utrecht Publishing \& Archiving Service},
+ author = {Oerlemans, Johannes},
+ year = {2008},
+ note = {OCLC: 600830356}
+}
+
+@article{tarasov_data-calibrated_2012,
+ series = {Sea {Level} and {Ice} {Sheet} {Evolution}: {A} {PALSEA} {Special} {Edition}},
+ title = {A data-calibrated distribution of deglacial chronologies for the {North} {American} ice complex from glaciological modeling},
+ volume = {315-316},
+ issn = {0012-821X},
+ url = {http://www.sciencedirect.com/science/article/pii/S0012821X11005243},
+ doi = {10.1016/j.epsl.2011.09.010},
+ abstract = {Past deglacial ice sheet reconstructions have generally relied upon discipline-specific constraints with no attention given to the determination of objective confidence intervals. Reconstructions based on geophysical inversion of relative sea level (RSL) data have the advantage of large sets of proxy data but lack ice-mechanical constraints. Conversely, reconstructions based on dynamical ice sheet models are glaciologically self-consistent, but depend on poorly constrained climate forcings and sub-glacial processes. As an example of a much better constrained methodology that computes explicit error bars, we present a distribution of high-resolution glaciologically-self-consistent deglacial histories for the North American ice complex calibrated against a large set of RSL, marine limit, and geodetic data. The history is derived from ensemble-based analyses using the 3D MUN glacial systems model and a high-resolution ice-margin chronology derived from geological and geomorphological observations. Isostatic response is computed with the VM5a viscosity structure. Bayesian calibration of the model is carried out using Markov Chain Monte Carlo methods in combination with artificial neural networks trained to the model results. The calibration provides a posterior distribution for model parameters (and thereby modeled glacial histories) given the observational data sets that takes data uncertainty into account. Final ensemble results also account for fits between computed and observed strandlines and marine limits. Given the model (including choice of calibration parameters), input and constraint data sets, and VM5a earth rheology, we find the North American contribution to mwp1a was likely between 9.4 and 13.2m eustatic over a 500year interval. This is more than half of the total 16 to 26m meltwater pulse over 500 to 700years (with lower values being more probable) indicated by the Barbados coral record (Fairbanks, 1989, Peltier and Fairbanks, 2006) if one assumes a 5meter living range for the Acropora Palmata coral. 20ka ice volume for North America was likely 70.1±2.0m eustatic, or about 60\% of the total contribution to eustatic sea level change. We suspect that the potentially most critical unquantified uncertainties in our analyses are those related to model structure (especially climate forcing), deglacial ice margin chronology, and earth rheology.},
+ language = {en},
+ urldate = {2020-04-28},
+ journal = {Earth and Planetary Science Letters},
+ author = {Tarasov, Lev and Dyke, Arthur S. and Neal, Radford M. and Peltier, W. R.},
+ month = jan,
+ year = {2012},
+ keywords = {uncertainty, glacial model, ice sheet reconstruction, Laurentide deglaciation, meltwater pulse, model calibration},
+ pages = {30--40},
+ file = {Tarasov et al_2012_A data-calibrated distribution of deglacial chronologies for the North American.pdf:/home/durbank/ODrive/Zotero attachments/Tarasov et al_2012_A data-calibrated distribution of deglacial chronologies for the North American.pdf:application/pdf}
+}
+
+@article{zemp_distributed_2007,
+ series = {Climate {Change} {Impacts} on {Mountain} {Glaciers} and {Permafrost}},
+ title = {Distributed modelling of the regional climatic equilibrium line altitude of glaciers in the {European} {Alps}},
+ volume = {56},
+ issn = {0921-8181},
+ url = {http://www.sciencedirect.com/science/article/pii/S0921818106001585},
+ doi = {10.1016/j.gloplacha.2006.07.002},
+ abstract = {Glaciers are among the key indicators of ongoing climate change. The equilibrium line altitude is a theoretical line which defines the altitude at which annual accumulation equals the ablation. It represents the lowest boundary of the climatic glacierisation and, therefore, is an excellent proxy for climate variability. In this study we introduce a simple approach for modelling the glacier distribution at high spatial resolution over entire mountain ranges using a minimum of input data. An empirical relationship between precipitation and temperature at the steady-state equilibrium line altitude (ELA0), is derived from direct glaciological mass balance measurements. Using geographical information systems (GIS) and a digital elevation model, this relationship is then applied over a spatial domain, to a so-called distributed modelling of the regional climatic ELA0 (rcELA0) and the climatic accumulation area (cAA) of 1971–1990 over the entire European Alps. A sensitivity study shows that a change in rcELA0 of ±100 m is caused by a temperature change of ±1 °C or a precipitation decrease of 20\% and increase of 27\%, respectively. The modelled cAA of 1971–1990 agrees well with glacier outlines from the 1973 Swiss Glacier Inventory. Assuming a warming of 0.6 °C between 1850 and 1971–1990 leads to a mean rcELA0 rise of 75 m and a corresponding cAA reduction of 26\%. A further rise in temperature of 3 °C accompanied by an increase in precipitation of 10\% leads to a further mean rise of the rcELA0 of about 340 m and reduces the cAA of 1971–1990 by 74\%.},
+ language = {en},
+ number = {1},
+ urldate = {2020-04-28},
+ journal = {Global and Planetary Change},
+ author = {Zemp, Michael and Hoelzle, Martin and Haeberli, Wilfried},
+ month = mar,
+ year = {2007},
+ keywords = {climate change, climate at equilibrium line altitude, geographical information systems, glacier},
+ pages = {83--100},
+ file = {Zemp et al_2007_Distributed modelling of the regional climatic equilibrium line altitude of.pdf:/home/durbank/ODrive/Zotero attachments/Zemp et al_2007_Distributed modelling of the regional climatic equilibrium line altitude of.pdf:application/pdf}
+}
+
+@book{cuffey_physics_2010,
+ title = {The {Physics} of {Glaciers}},
+ isbn = {978-0-08-091912-6},
+ abstract = {The Physics of Glaciers, Fourth Edition, discusses the physical principles that underlie the behavior and characteristics of glaciers. The term glacier refers to all bodies of ice created by the accumulation of snowfall, e.g., mountain glaciers, ice caps, continental ice sheets, and ice shelves. Glaciology—the study of all forms of ice—is an interdisciplinary field encompassing physics, geology, atmospheric science, mathematics, and others. This book covers various aspects of glacier studies, including the transformation of snow to ice, grain-scale structures and ice deformation, mass exchange processes, glacial hydrology, glacier flow, and the impact of climate change. The present edition features two new chapters: “Ice Sheets and the Earth System and “Ice, Sea Level, and Contemporary Climate Change. The chapter on ice core studies has been updated from the previous version with new material. The materials on the flow of mountain glaciers, ice sheets, ice streams, and ice shelves have been combined into a single chapter entitled “The Flow of Ice Masses. Completely updated and revised, with 30\% new material including climate changeAccessible to students, and an essential guide for researchersAuthored by preeminent glaciologists},
+ language = {en},
+ publisher = {Academic Press},
+ author = {Cuffey, Kurt M. and Paterson, W. S. B.},
+ month = jun,
+ year = {2010},
+ note = {Google-Books-ID: Jca2v1u1EKEC},
+ keywords = {Science / Earth Sciences / Geology, Science / Physics / Geophysics}
+}
+
+@phdthesis{farinotti_simple_2010,
+ type = {Doctoral {Thesis}},
+ title = {Simple methods for inferring glacier ice-thickness and snow-accumulation distribution},
+ copyright = {http://rightsstatements.org/page/InC-NC/1.0/},
+ url = {https://www.research-collection.ethz.ch/handle/20.500.11850/152254},
+ language = {en},
+ urldate = {2020-04-29},
+ school = {ETH Zurich},
+ author = {Farinotti, Daniel},
+ year = {2010},
+ doi = {10.3929/ethz-a-006194634},
+ note = {Accepted: 2017-06-13T11:46:19Z},
+ file = {Full Text PDF:/home/durbank/Zotero/storage/CNJD25WI/Farinotti - 2010 - Simple methods for inferring glacier ice-thickness.pdf:application/pdf}
+}
+
+@techreport{wgms_fluctuations_2019,
+ address = {Zurich, Switzerland},
+ title = {Fluctuations of {Glaciers} {Database}},
+ url = {http://dx.doi.org/10.5904/wgms-fog-2019-12},
+ institution = {World Glacier Monitoring Service (WGMS)},
+ author = {WGMS},
+ year = {2019}
+}
+
+@mastersthesis{keeler_development_2015,
+ address = {Provo, Utah USA},
+ title = {Development and {Validation} of a {Physically} {Based} {ELA} {Model} and its {Application} to the {Younger} {Dryas} {Event} in the {Graubünden} {Alps}, {Switzerland}},
+ url = {https://scholarsarchive.byu.edu/etd/6142},
+ abstract = {The rapid rate of global warming currently underway highlights the need for a deeper understanding of abrupt climate change. The Younger Dryas is a Late-Glacial climate event of widespread and unusually rapid change whose study can help us address this need for increased understanding. Reconstructions from the glacial record offer important contributions to our understanding of the Younger Dryas due to (among other things) the direct physical response of glaciers to even minor perturbations in climate. Because the glacier equilibrium line altitude (ELA) provides a more explicit comparison of climate than properties such as glacier length or area, ELA methods lend themselves well to paleoclimate applications and allow for more direct comparisons in space and time. Here we present a physically based ELA model for alpine paleoglacier climate reconstructions that accounts for differences in glacier width, glacier shape, bed topography and ice thickness, and includes error estimates using Monte Carlo simulations. We validate the ELA model with published mass balance measurements from 4 modern glaciers in the Swiss Alps. We then use the ELA model, combined with a temperature index model, to estimate the changes in temperature and precipitation between the Younger Dryas (constrained by 10Be surface exposure ages) and the present day for three glacier systems in the Graubünden Alps. Our results indicate an ELA depression in this area of 320 m ±51 m during the Younger Dryas relative to today. This ELA depression represents annual mean temperatures 2.29 °C ±1.32 °C cooler relative to today in the region, which corresponds to a decrease in mean summer temperatures of 1.47 °C ±0.73 °C. Our results indicate relatively small changes in summer temperature dominate over other climate changes for the Younger Dryas paleoglaciers in the Alps. This ELA-based paleoclimate reconstruction offers a simple, fast, and cost-effective alternative to many other paleoclimate reconstruction methods. Continued application of the ELA model to more regions will lead to an improved understanding of the Younger Dryas in the Alps, and by extension, of rapid climate events generally.},
+ school = {Brigham Young University},
+ author = {Keeler, Durban},
+ month = nov,
+ year = {2015}
+}
+
+@article{thomas_sensitivity_2014,
+ title = {Sensitivity of digital elevation models: {The} scenario from two tropical mountain river basins of the {Western} {Ghats}, {India}},
+ volume = {5},
+ issn = {1674-9871},
+ shorttitle = {Sensitivity of digital elevation models},
+ url = {http://www.sciencedirect.com/science/article/pii/S1674987114000036},
+ doi = {10.1016/j.gsf.2013.12.008},
+ abstract = {The paper evaluates sensitivity of various spaceborne digital elevation models (DEMs), viz., Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Shuttle Radar Topography Mapping Mission (SRTM) and Global Multi-resolution Terrain Elevation Data 2010 (GMTED), in comparison with the DEM (TOPO) derived from contour data of 20 m interval of Survey of India topographic sheets of 1: 50,000 scale. Several topographic attributes, such as elevation (above mean sea level), relative relief, slope, aspect, curvature, slope-length and -steepness (LS) factor, terrain ruggedness index (TRI), topographic wetness index (TWI), hypsometric integral (Ihyp) and drainage network attributes (stream number and stream length) of two tropical mountain river basins, viz., Muthirapuzha River Basin and Pambar River Basin are compared to evaluate the variations. Though the basins are comparable in extent, they differ in respect of terrain characteristics and climate. The results suggest that ASTER and SRTM provide equally reliable representation of topography portrayed by TOPO and the topographic attributes extracted from the spaceborne DEMs are in agreement with those derived from TOPO. Despite the coarser resolution, SRTM shows relatively higher vertical accuracy (RMSE = 23 and 20 m respectively in MRB and PRB) compared to ASTER (RMSE = 33 and 24 m) and GMTED (RMSE = 59 and 48 m). Vertical accuracy of all the spaceborne DEMs is influenced by relief of the terrain as well as type of vegetation. Further, GMTED shows significant deviation for most of the attributes, indicating its inability for mountain-river-basin-scale studies.},
+ language = {en},
+ number = {6},
+ urldate = {2020-05-04},
+ journal = {Geoscience Frontiers},
+ author = {Thomas, Jobin and Joseph, Sabu and Thrivikramji, K. P. and Arunkumar, K. S.},
+ month = nov,
+ year = {2014},
+ keywords = {ASTER, DEM, GMTED, SRTM, Tropical mountain river basins, Western Ghats},
+ pages = {893--909},
+ file = {Thomas et al_2014_Sensitivity of digital elevation models.pdf:/home/durbank/ODrive/Zotero attachments/Thomas et al_2014_Sensitivity of digital elevation models.pdf:application/pdf}
+}
+
+@techreport{aster_aster_2009,
+ title = {{ASTER} {Global} {DEM} {Validation} {Summary} {Report}},
+ url = {https://lpdaac.usgs.gov/documents/28/ASTER_GDEM_Validation_1_Summary_Report.pdf},
+ abstract = {This report presents highlights and the most relevant results from initial studies conducted to validate and otherwise characterize the new global digital elevation model (DEM) produced by the Ministry of Economy, Trade and Industry (METI) of Japan and the United States National Aeronautics and Space Administration (NASA) from optical stereo data acquired by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER). The ASTER Global DEM (GDEM) was released to the public on June 29, 2009. Results and findings presented in this report were instrumental in the decision by METI and NASA to release the ASTER GDEM.},
+ institution = {ASTER GDEM Validation Team: METI, NASA and USGS in cooperation with NGA and other collaborators},
+ author = {ASTER},
+ year = {2009},
+ pages = {28},
+ file = {DTED_2009_ASTER Global DEM Validation Summary Report.pdf:/home/durbank/ODrive/Zotero attachments/DTED_2009_ASTER Global DEM Validation Summary Report.pdf:application/pdf}
+}
+
+
+@article{braithwaite_estimating_2009,
+ title = {Estimating equilibrium-line altitude ({ELA}) from glacier inventory data},
+ volume = {50},
+ issn = {0260-3055, 1727-5644},
+ url = {https://www.cambridge.org/core/journals/annals-of-glaciology/article/estimating-equilibriumline-altitude-ela-from-glacier-inventory-data/89DEF6B31BA17C0DF6B5B566BCD2CC64},
+ doi = {10.3189/172756410790595930},
+ abstract = {A glacier’s most fundamental altitude is the equilibrium-line altitude (ELA) because it divides the glacier into ablation and accumulation areas. The best parameterization of the ELA for glacier inventory is the balanced-budget ELA. We discuss direct estimation of balanced-budget ELA from mass-balance data for individual glaciers, and indirect estimation of balanced-budget ELA from simple topographic parameters available from the World Glacier Inventory (WGI), i.e. the area-median and maximum and minimum altitudes. Mass balance and ELA for individual glaciers are usually strongly correlated and we calculate balanced-budget ELA from the regression equation linking the two. We then compare balanced-budget ELA with area-median and mid-range altitudes for the 94 glaciers for which we have all the necessary data. The different ELA estimates agree well enough (±82 to ±125 m) to describe geographical variations in ELA and for application of glacier–climate models to glacier inventory data. Mid-range and area-median altitudes are already available for tens of thousands of glaciers in the current WGI and should be evaluated in future inventories.},
+ language = {en},
+ number = {53},
+ urldate = {2020-11-11},
+ journal = {Annals of Glaciology},
+ author = {Braithwaite, R. J. and Raper, S. C. B.},
+ year = {2009},
+ note = {Publisher: Cambridge University Press},
+ pages = {127--132}
+}
+
+@article{benn_exceltm_2010,
+ title = {An {ExcelTM} spreadsheet program for reconstructing the surface profile of former mountain glaciers and ice caps},
+ volume = {36},
+ issn = {0098-3004},
+ url = {http://www.sciencedirect.com/science/article/pii/S0098300410000531},
+ doi = {10.1016/j.cageo.2009.09.016},
+ abstract = {A new Excel™ spreadsheet program is introduced, which calculates the surface profiles of former glaciers using an exact solution of a ‘perfectly plastic’ glacier model. Two versions of the model are presented. The basic model requires only bed topography along a flowline and a yield stress for ice as inputs. The latter can be tuned using ‘target ice elevations’ derived from geomorphological mapping. A more sophisticated form of the model allows the yield stress to vary along the flowline, and incorporates the effect of valley-side drag on the glacier profile.},
+ language = {en},
+ number = {5},
+ urldate = {2020-11-11},
+ journal = {Computers \& Geosciences},
+ author = {Benn, Douglas I. and Hulton, Nicholas R. J.},
+ month = may,
+ year = {2010},
+ keywords = {Glacier reconstruction, Ice sheet surface, Shape factor},
+ pages = {605--610}
+}
+
+@article{pellitero_gis_2015,
+ title = {A {GIS} tool for automatic calculation of glacier equilibrium-line altitudes},
+ volume = {82},
+ issn = {0098-3004},
+ url = {http://www.sciencedirect.com/science/article/pii/S0098300415001090},
+ doi = {10.1016/j.cageo.2015.05.005},
+ abstract = {A toolbox for the automated calculation of glacier equilibrium-line altitudes (ELAs) using the Accumulation Area Ratio, Area-Altitude Balance Ratio, Area-Altitude and Kurowski methods is presented. These are the most commonly-used methods of ELA calculation in palaeo-glacier reconstructions. The toolbox has been coded in Python and runs in ArcGIS requiring only the reconstructed surface of the palaeo-glacier (a DEM) as input. Through fast and automatic calculation this toolbox simplifies the process of ELA determination and can successfully work both for a single glacier and for large datasets of multiple glaciers.},
+ language = {en},
+ urldate = {2020-11-11},
+ journal = {Computers \& Geosciences},
+ author = {Pellitero, Ramón and Rea, Brice R. and Spagnolo, Matteo and Bakke, Jostein and Hughes, Philip and Ivy-Ochs, Susan and Lukas, Sven and Ribolini, Adriano},
+ month = sep,
+ year = {2015},
+ keywords = {GIS, Accumulation area balance ratio, Accumulation area ratio, Equilibrium-line altitude, Glacier, Python},
+ pages = {55--62}
+}
+
+@article{pellitero_glare_2016,
+ title = {{GlaRe}, a {GIS} tool to reconstruct the {3D} surface of palaeoglaciers},
+ volume = {94},
+ issn = {0098-3004},
+ url = {http://www.sciencedirect.com/science/article/pii/S0098300416301558},
+ doi = {10.1016/j.cageo.2016.06.008},
+ abstract = {Glacier reconstructions are widely used in palaeoclimatic studies and this paper presents a new semi-automated method for generating glacier reconstructions: GlaRe, is a toolbox coded in Python and operating in ArcGIS. This toolbox provides tools to generate the ice thickness from the bed topography along a palaeoglacier flowline applying the standard flow law for ice, and generates the 3D surface of the palaeoglacier using multiple interpolation methods. The toolbox performance has been evaluated using two extant glaciers, an icefield and a cirque/valley glacier from which the subglacial topography is known, using the basic reconstruction routine in GlaRe. Results in terms of ice surface, ice extent and equilibrium line altitude show excellent agreement that confirms the robustness of this procedure in the reconstruction of palaeoglaciers from glacial landforms such as frontal moraines.},
+ language = {en},
+ urldate = {2020-11-11},
+ journal = {Computers \& Geosciences},
+ author = {Pellitero, Ramón and Rea, Brice R. and Spagnolo, Matteo and Bakke, Jostein and Ivy-Ochs, Susan and Frew, Craig R. and Hughes, Philip and Ribolini, Adriano and Lukas, Sven and Renssen, Hans},
+ month = sep,
+ year = {2016},
+ keywords = {GIS, Glacier reconstruction, Python, Flowline equilibrium profile},
+ pages = {77--85}
+}
diff --git a/Docs/ELA-model_tex/main.tex b/Docs/ELA-model_tex/main.tex
new file mode 100644
index 0000000..1f3a3a7
--- /dev/null
+++ b/Docs/ELA-model_tex/main.tex
@@ -0,0 +1,418 @@
+\listfiles
+\documentclass[review]{elsarticle}
+
+\usepackage{lineno,hyperref}
+\usepackage{amsmath}
+\modulolinenumbers[5]
+\usepackage{hyperref}
+\def\UrlBreaks{\do\/\do-}
+\usepackage[margin=1in]{geometry}
+
+\newcommand{\mytilde}{\raise.17ex\hbox{$\scriptstyle\mathtt{\sim}$}}
+
+\journal{MethodsX}
+
+\begin{document}
+
+\begin{frontmatter}
+
+\title{A first-order flexible ELA model based on geomorphic constraints\tnoteref{mytitlenote}}
+\tnotetext[mytitlenote]{Co-submission to MethodsX with JQSR\_2020\_146}
+
+\author[Utah]{Durban G. Keeler\corref{corresponding}}
+\cortext[mycorrespondingauthor]{Durban Keeler}
+\ead{durban.keeler@gmail.com}
+\author[Utah]{Summer Rupper}
+\address[Utah]{Department of Geography, University of Utah, Salt Lake City, UT USA}
+
+\author[Columbia]{Joerg Schaefer}
+\address[Columbia]{Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY USA}
+
+\begin{abstract}
+Alpine glaciers, with their valuable combination of highly sensitive response to climate and near-global extent, are powerful tools for investigating previous and present climate changes.
+They also represent critical water resources for areas around the globe, with the potential for far-reaching effects in a warming world.
+Advancements to understand and model glacial changes and the variables influencing them are therefore paramount.
+Many glacier models fall into one of two endmembers: either highly complex transient models requiring careful tuning of multiple parameters to individual glaciers, or basic empirical correlations of glacier area and length with few considerations for local and regional variations in characteristics.
+Here we detail a physical steady-state model for alpine glaciers relating directly to glacier mass balance (via the equilibrium line altitude) while retaining the simplicity of other morphology methods, and simultaneously including error estimates.
+We provide custom MATLAB functions as a user-friendly and generally-applicable method to estimate glacier equilibrium line altitudes from only a limited number of glacier bed topography and glacier width measurements.
+As a test of the model’s efficacy, we compare the model results for present-day glaciers in the Swiss Alps with previously published estimates of equilibrium line altitudes and intermediate model outputs.
+
+\begin{itemize}
+ \item The method estimates glacier equilibrium line altitudes from a limited set of bed topography measurements and constraints on glacier width
+ \item The method is based on continuity equations, reducing the need for empirical coefficients tuned with measured data
+ \item The method uses Monte Carlo sampling and bootstrapping to generate uncertainty bounds on the equilibrium line altitude estimates
+\end{itemize}
+
+\end{abstract}
+
+\begin{keyword}
+glacier modeling\sep climate reconstruction\sep equilibrium line altitude
+\end{keyword}
+
+\begin{figure}
+ \centering
+ \includegraphics[width=\textwidth]{Figures/graph-abstract.png}
+ \caption{Graphical abstract describing the ELA modeling process used in this study.}
+ \label{fig:abstract}
+\end{figure}
+
+\end{frontmatter}
+
+\linenumbers
+
+% \section*{Graphical Abstract}
+
+\section{Introduction}\label{sec:intro}
+
+Changes in glacier size and extent are fundamentally related to the mass balance of the glacier.
+An annual mass surplus (when net accumulation exceeds ablation) leads to glacier growth, while a deficit leads to glacier retreat.
+Such transitions can be visually drastic, with some glaciers changing by kilometers in response to minor perturbations in climate.
+A more comparable measure of the glacier response to climate change than glacier area or length changes is the concept of the equilibrium line altitude (ELA).
+The ELA is the boundary between the accumulation and ablation zones on a glacier and represents the elevation at which the annual amount of mass added through accumulation exactly equals the annual amount of mass lost through ablation \cite{cuffey_physics_2010}.
+As a direct measure of glacier mass balance, the ELA facilitates explicit comparisons of climate in space and time by accounting for dependencies on glacier size, extent, and shape, and by integrating the myriad variables that can drive changes in climate into a single comparable metric.
+
+The most accurate method to determine the ELA is direct measurements of the distribution of accumulation and ablation on a glacier.
+Such measurements, however, are only available for a small subset of modern glaciers, and entirely absent for paleoglaciers.
+Several empirical proxy methods have therefore been developed to estimate ELA when direct measurements of mass balance are not possible.
+Some of the more widely used are the accumulation area ratio (AAR), the toe-to-headwall altitude ratio (THAR), the area-weighted mean altitude (AMA), and the median elevation of the glacier (MEG) \cite{benn_mass_2000, braithwaite_estimating_2009}.
+Such statistical methods are useful within many contexts.
+\cite{braithwaite_estimating_2009} for instance found an overall $R^2 = 0.99$ for 94 glaciers between steady-state ELA and one of the simplest ways to estimate ELA (MEG).
+Estimates of individual glaciers, however, can significantly diverge.
+Moreover, the empirical nature of these methods make them relatively easy to apply, but is also one of the primary inherent limitations on their applicability and interpretability.
+Because the empirical coefficients are derived from aggregate glacier data sets, these models are only valid for glaciers within the boundary conditions of the training data, typically with no a priori method to determine whether such an assumption is valid when applied to other regions or other periods of time.
+Such regionally-tuned coefficients can vary widely.
+AAR values among modern glaciers, for instance, vary between 0.22 and 0.72, resulting in a large range in possible ELA estimates depending on the value chosen \cite{zemp_distributed_2007}.
+The other empirical methods for ELA reconstructions suffer from similar concerns as the AAR method.
+Additionally, such methods provide no insight into the sensitivity or uncertainty of the estimates to glacier characteristics such as bed topography or areal distribution.
+All of these methods, although useful in many circumstances, highlight the need for additional progress to help better constrain ELA estimates in a robust, self-consistent manner and place the estimate within the context of model error, while still requiring minimum inputs.
+
+This manuscript presents a method to reconstruct ELA estimates based on continuity equations, while only requiring estimates of bed topography, glacier length, and glacier width.
+The ELA model also generates intermediate results of continuous modeled bed topography, ice surface elevation, and glacier width along the length of the glacier (Figure 1).
+An added strength of this model is that such intermediate outputs allow for increased diagnostics on model performance or troubleshooting.
+The required inputs for the model are similar to those required for simple geomorphic methods, but based on continuity equations rather than relying purely on empirical correlations while also accounting for physical errors in measurements.
+This method can equally apply to existing glaciers or paleo-glacier extents where glacial moraines are adequately preserved.
+
+\section{Methodology}\label{sec:methods}
+
+The presented model is largely derived from a simple linear glacier-length model presented in \cite{oerlemans_minimal_2008}, with modifications specific to quantifying ELAs and ELA changes.
+The limited model inputs necessarily require simplifying assumptions that do not include all details pertinent to specific glaciers.
+Such details can be significant for some applications (e.g. dynamic modeling of glacier response, higher order surface energy and mass balance modeling, etc.), and other methods would be better suited to these circumstances.
+The proposed model is specifically intended to estimate the ELA of snow-fed, clean ice, temperate valley glaciers with relatively simple bed and areal geometries.
+
+The ELA model also provides analytical constraints on the error associated with model outputs.
+Such uncertainties are fundamental in determining the significance and reliability of results, but rigorous physical uncertainties of ELA estimates are rarely presented in paleo-glacier research, either because such uncertainties are difficult to assign for geomorphic methods like THAR and AAR or because higher order models are sufficiently complex to challenge error propagation.
+Uncertainty estimates in this study are calculated based on Monte Carlo simulations of bootstrapped residuals of the input parameters.
+These uncertainties give insight into the range of plausible ELA values based on both uncertainty of input parameters and the ability of the model assumptions to accurately represent those inputs.
+
+Other computational methods exist to estimate the ELA of paleoglaciers.
+\cite{benn_exceltm_2010} presents an Excel\textsuperscript{TM} spreadsheet to calculate the ice surface profiles of a former mountain glacier or ice cap, given bed topography and a yield stress.
+\cite{pellitero_gis_2015} provides a Python-based ArcGIS toolbox to automatically calculate glacier ELAs with a choice of methods (Accumulation-Area Ratio, Area-Altitude-Balance Ratio, Area-Altitude, or median elevation).
+The tool requires a DEM of the reconstructed glacier surface as input.
+\cite{pellitero_glare_2016} builds on \cite{pellitero_gis_2015} by adding the ability to estimate and interpolate a paleoglacier ice surface given the 3D bed topography and a center flowline.
+These methods represent important steps forward by incorporating ice flow laws and automating much of the process in ELA calculations in an accessible manner.
+
+Similar to \cite{benn_exceltm_2010, pellitero_glare_2016}, the presented ELA model also estimates the glacier surface based on centerline ice flow, given bed topography.
+It also aims to automate many of the steps in calculating an ELA to provide an easy-to-use and widely applicable method of ELA estimation.
+The main advantages of the presented ELA model are (1) less reliance on empirical relations which require tuning of coefficients and (2) the inclusion of robust uncertainty bounds based on Monte Carlo sampling.
+The decreased reliance on empirical correlations permits more widespread applicability without a priori knowledge and tuning of model coefficients.
+Such coefficients are determined regionally using glacier mass balance measurements from many glaciers in the area, even though nearby glacier values can diverge notably \cite{benn_mass_2000}.
+The proposed ELA model seeks to avoid these empirical complications by estimating the ELA directly from continuity equations.
+Another important advantage of the proposed ELA model is the incorporation of Monte Carlo bootstrapping.
+This permits estimation of robust uncertainties both from error in input parameters and from model limitations.
+These uncertainties give proper context to how accurate results are for specific glaciers, as well as how sensitive these glaciers are to various inputs.
+Although this requires more complex modeling and computation to generate results compared to empirical methods, most of this complexity is fully automated and abstracted away from the user.
+It therefore retains simplicity, applicability, and usability while implementing more complex modeling and more robust uncertainty estimates.
+
+\subsection{Balance equation}
+
+The fundamental basis of the ELA model is an integrated balance formula (Equation \ref{eq:balance}) for steady-state glaciers from \cite{oerlemans_minimal_2008},
+
+\begin{equation}
+\label{eq:balance}
+B_n = \int_0^L \dot{b}dx = \beta \int_0^L \left[ w(x) \left( H(x)+z(x)-ELA \right) \right]dx
+\end{equation}
+
+where $B_n$ is the total net balance, $x$ is the distance down glacier, $\dot{b}$ is the specific balance rate at $x$, $L$ is the glacier length, $\beta$ is the balance gradient, $w(x)$ is the glacier width at $x$, $H(x)$ the ice thickness at $x$, $z(x)$ represents the valley topography, and $ELA$ is the equilibrium line altitude.
+In steady state conditions (like we assume for glaciers with well-developed moraine sequences), the total net balance is zero.
+The balance gradient $\beta$ can be dropped in this case, and Equation \ref{eq:balance} can then be adapted to solve for the $ELA$ (Equation \ref{eq:ela}).
+
+\begin{equation}
+\label{eq:ela}
+ELA = \frac{\int_0^L w(x) H(x)dx + \int_0^L w(x) z(x)dx}{\int_0^L w(x) dx}
+\end{equation}
+
+We then estimate each of the three components ($H(x)$, $w(x)$, and $z(x)$) along the length of the glacier and solve for each component using trapezoidal numerical integration to derive an estimate for $ELA$.
+Methods for the estimation of each of these components are detailed below.
+
+\subsection{Glacier bed modeling}
+
+Bed topography measurements follow a representative 1D line along the glacier profile taken down the center of the glacier.
+We then estimate $z(x)$ from a best-fit two-term exponential curve of this 1D profile line (Equation \ref{eq:bed}), where $a$, $b$, $c$, and $d$ are fitting coefficients optimized in the model using the elevation data inputs (see Figure \ref{fig:hyp} for examples).
+Optimizations in this ELA model use nonlinear least squares regression based on \textit{trust region} algorithms \cite{more_computing_1983}.
+
+\begin{equation}
+\label{eq:bed}
+z(x) = ae^{bx} + ce^{dx}
+\end{equation}
+
+This two-term exponential estimate is best suited for valleys with relatively simple bed topographies.
+Caution should be used when applying this method to glacier beds with more complex bed features, such as steep cliffs or over deepenings, as these are not always readily captured in the model.
+
+\subsection{Ice thickness modeling}
+
+To first order, the thickness of a glacier depends largely on the slope and shear stress at the bed \cite{cuffey_physics_2010}.
+The simplest equation to approximate ice thickness is therefore
+
+\begin{equation}
+\label{eq:ice0}
+H = \frac{\tau}{\rho g \sin\theta}
+\end{equation}
+
+where $H$ is the ice thickness (m), $\tau$ is the basal shear stress (Pa), $\rho$ is the ice density ($\:kg/m^3$), $g$ is acceleration due to gravity ($\:m/s^2$), and $\theta$ is the angle at the bed interface \cite{cuffey_physics_2010}.
+In areas with shallow slopes, however, Equation \ref{eq:ice0} leads to ice thickness unrealistically approaching infinity.
+\cite{oerlemans_minimal_2008} demonstrates a square root relation between length and ice thickness (assuming perfect plasticity), which we incorporate into our estimates in order to address this issue.
+
+\begin{equation}
+\label{eq:ice-thick}
+H_m = \frac{2}{3} \sqrt{\frac{\tau L}{\rho g \left( 1+\sin\theta \right)}}
+\end{equation}
+
+Equation \ref{eq:ice-thick}, however, gives the mean ice thickness ($H_m$) for the glacier, rather than continuous values along its length.
+To model ice thickness profiles, we assume a parabolic distribution (true of a perfectly plastic glacier on a flat bed) around the mean ice thickness (see Figure \ref{fig:ice} for examples).
+The basal shear stress ($\tau$) is assumed to scale with ice thickness, following the relationship presented in \cite{haeberli_application_1995} (Equation \ref{eq:tau}), where $\Delta z$ is the difference between the minimum and maximum bed elevation.
+
+\begin{equation}
+\label{eq:tau}
+\begin{split}
+\Delta z > 1600 \:m \Longrightarrow \tau = 150 \:kPa \\
+500 \:m \le \Delta z \le 1600 \:m \Longrightarrow \tau = 0.005 + 1.598\Delta z - 0.435\Delta z^2 \:kPa \\
+\Delta z < 500 \:m \Longrightarrow \tau = 3\Delta z \:kPa
+\end{split}
+\end{equation}
+
+\subsection{Glacier width modeling}
+
+Due to the high diversity in glacier shape/geometry, estimating the plan-view profile of the glacier in a consistent yet simple manner is difficult.
+Additionally, accurately constraining the width of the accumulation area for paleoglaciers presents further challenges, due to a lack of preserved moraines or other features delineating glacier boundaries in these areas.
+To best cope with these issues, we estimate glacier width using an exponential formula (Equation \ref{eq:width}) of the same form as presented in \cite{oerlemans_minimal_2008}.
+The initial starting parameters are the minimum width of the glacier at the toe ($w0$), maximum glacier width in the accumulation zone ($w_{max}$), the distance down glacier ($x$), and the distance down glacier to the point of maximum width ($L_{Wmax}$).
+
+\begin{equation}
+\label{eq:width}
+w(x) = w_0 + \frac{w_{max}-w_0}{L_{Wmax}}xe^{1-\frac{x}{L_{Wmax}}}
+\end{equation}
+
+This produces an exponential curve, following the general shape of many glaciers.
+The model then performs least squares nonlinear curve fitting (again based on trust region techniques) on the initial parameter estimates to optimize the fit to the input width estimates (see Figure \ref{fig:width} for examples).
+
+The model can also incorporate glacier tributaries.
+The tributaries are initially modeled as independent glaciers, including profile centerline elevations and width measurements.
+The calculated tributary glacier volume is then added to the main glacier at corresponding elevation levels as additional modeled glacier width.
+Added caution should be exercised with this model when including tributary glaciers, as the glacier plan profile can depart more severely from the assumed idealized shape constraints.
+
+\subsection{Monte Carlo simulations}
+
+We perform Monte Carlo simulations to capture the distribution of plausible ELAs for a given glacier.
+Such estimation of uncertainty is important to adequately compare the significance of results, particularly if attempting to compare results from differing methodologies or between regions.
+Monte Carlo methods are widely used in modeling of glacier mass and energy balance for uncertainty estimation \cite{machguth_exploring_2008}.
+In our approach, we perform bootstrapping with replacement techniques to incorporate the uncertainty of input parameters and to include any known errors in those parameters, assuming Gaussian error distributions.
+Each model run consists of 1,000 simulations in order to approximate a continuous distribution in plausible ELA values.
+The Monte Carlo simulations do increase the computational load, especially compared to the automated methods of \cite{benn_exceltm_2010, pellitero_gis_2015}, taking \mytilde1 minute to process one glacier on a single core.
+The model code, however, utilizes parallel processing, enabling much greater scalability to larger data sets with the proper hardware.
+
+\section{Data and analysis workflow}
+
+The complete ELA model MATLAB code is publicly available (\url{https://github.com/durbank/ELA-model}), with v0.1.0 the particular version used in this manuscript.
+% Detailed documentation on using the ELA model can also be found in the GitHub repository, as well as an example script demonstrating the model on four glacier test sites (see the Model Validation section for details).
+In brief, the ELA model function ELA\_calc.m requires two dataset inputs (discrete estimates of bed topography and discrete estimates of glacier width, both measured downglacier along the centerline of the glacier valley) and the number of Monte Carlo simulations to perform.
+Approximately ten quasi-equally spaced points along the length of the glacier are often sufficient, though the optimum number depends on the length and complexity of the bed topography and glacier geometry.
+To avoid issues of model extrapolation (and to automatically include the overall glacier length), both the toe and the head of the glacier should be included in these measurements.
+The ELA model input data should be provided as a MATLAB structure with four fields, as summarized in Table \ref{tab:inputs}.
+Tributary glaciers, if present, should be provided as variable input arguments (formatted as MATLAB structures according to Table \ref{tab:inputs}) after the number of simulations to perform.
+The format\_inputs.m function takes .csv files of glacier bed topography and glacier width measurements and creates a properly-formatted MATLAB structure to serve as input to the ELA model.
+
+\begin{table}[]
+ \centering
+ \begin{tabular}{c|c|p{11cm}}
+ Field name & Dimensions & Field description \\
+ \hline
+ X\_dist & $[N \times 1]$ & Vector of glacier length from $0:N$, where N is the total length of the glacier in meters \\
+ Bed\_pts & $[n \times 2]$ & A matrix with positions along the glacier centerline (in meters) in the first column and corresponding bed elevation measurements (meters a.s.l.) in the second \\
+ Ice\_surf & $[n \times 2]$ & A matrix with positions along the glacier centerline (in meters) in the first column (this should match the first column in 'Bed\_pts') and corresponding ice surface elevation measurements (meters a.s.l.) in the second \\
+ Width\_pts & $[m \times 2]$ & A matrix with positions along the glacier centerline (in meters) in the first column and corresponding glacier width measurements (meters) in the second (widths should orthogonally intersect the centerline)
+ \end{tabular}
+ \caption{Required format for ELA model inputs}
+ \label{tab:inputs}
+\end{table}
+
+In addition to the inputs, there are model parameter assumptions built into the model prescribing the assumed errors for Monte Carlo sampling.
+Updating these assumptions to better reflect specific input data is a simple matter of editing the assigned values.
+Table \ref{tab:error} shows a summary of these parameters and their default values.
+
+\begin{table}[]
+ \centering
+ \begin{tabular}{c|c|p{10cm}}
+ Variable name & Default value & Variable description \\
+ \hline
+ zSTD & $25 \:meters$ & Standard deviation in measured glacier bed elevation \\
+ wSTD & $60 \:meters$ & Standard deviation in measured glacier width values \\
+ tau\_STD & $5.0 \times 10^4 \:Pa$ & Standard deviation in estimated basal shear stress (used in ice thickness calculations) \\
+ rho & $917 \:kg/m^3$ & Density of ice (used in ice thickness calculations) \\
+ g & $9.8 \:m/s^2$ & Acceleration due to gravity (used in ice thickness calculations)
+ \end{tabular}
+ \caption{ELA model error assumptions}
+ \label{tab:error}
+\end{table}
+
+For the development and validation of this model, we used a particular ArcGIS software workflow to generate the ELA model inputs.
+We include this workflow as a diagram (Figure \ref{fig:arc-flow}), but model inputs can be generated and provided using any desired methods, as long as they are properly formatted.
+
+\begin{figure}
+ \centering
+ \includegraphics[width=\textwidth]{Figures/arc-flow.png}
+ \caption{Flowchart showing the ArcGIS workflow used to generate ELA model inputs. Orange denotes steps performed in ArcMap, while blue denotes steps performed in MATLAB. The first step is to generate a characteristic centerline for the glacier. This centerline can be drawn freehand or calculated in some other way. Then extract the bed elevation along the centerline using the DEM input (recorded as distance along the centerline), and save to a temporary .xls file. Import this file into MATLAB and use the "ice\_thick.m" function (part of the ELA model) to estimate the ice surface elevation at points along the centerline transect. Import the ice surface results back to ArcMap and combine with the centerline transect values. The final steps require an outline of the glacier in question. These boundaries can be drawn from the aerial imagery for modern glaciers, or else from the moraines of paleoglaciers. In the case of paleoglacier moraines, the accumulation region of the glacier can be broadly defined by the valley boundaries. Calculate polylines at each discrete point along the transect at the elevation of the ice surface and orthogonal to the centerline at that point. The intersection of these orthogonal lines (at the elevation of the ice surface at that point) with either the glacier boundaries (in the case of modern glacier outlines and portions of paleoglaciers constrained by moraines) or the bed topography (in the case of the accumulation zone of paleoglaciers) defines the glacier width at each transect point. The results of distance down the glacier centerline, estimated bed elevation, and estimated glacier widths are then saved as a .csv file for import to the ELA model.}
+ \label{fig:arc-flow}
+\end{figure}
+
+\section{Model validation}
+
+We validate the ELA model by matching our reconstructions with direct observations of four modern glaciers in the European Alps.
+These glaciers were primarily selected due to the availability of data requisite for a data-model comparison (including present-day ice thickness, bed topography beneath the present-day glacier, mass balance measurements, aerial photography and DEMs).
+Moreover, the validation glaciers show variation in overall area, length, width, and elevation extent, thereby providing a range of possible glacier geometries (see Table \ref{tab:attrs} for differences).
+This range of glacier characteristics enables a more rigorous test of robustness and general applicability of the ELA model.
+Although all 4 validation glaciers are limited to the European Alps, the method should be more widely applicable to similar glaciers in other regions (clean-ice, land-terminating temperate mountain glaciers) without regional tuning of empirical coefficients.
+The four test glaciers are the Gries, Findel, Rhone, and Silvretta Glaciers.
+Three of these glaciers (Gries, Silvretta, and Findel) have continuous multi-year mass balance measurements from stake networks compiled by the World Glacier Monitoring Service (WGMS), and therefore make for the most compelling comparisons.
+The Rhone Glacier has mass balance measurements from a handful of isolated years, providing a less certain, but still useful comparison to the model and other glaciers.
+
+\begin{table}[]
+ \centering
+ \begin{tabular}{c|c|c|c|c|c|c|c}
+ Glacier & Min elevation & Max elevation & Observed ELA & Area & Length & Max width \\
+ \hline
+ Gries & $2500\:m$ & $3277\:m$ & $3084\:m$ & $4.28\:km^2$ & $5216\:m$ & $1441\:m$ \\
+ Silvretta & $2498\:m$ & $3078\:m$ & $2829\:m$ & $2.58\:km^2$ & $3193\:m$ & $1307\:m$ \\
+ Findel & $2600\:m$ & $3570\:m$ & $3233\:m$ & $14.5\:km^2$ & $7089\:m$ & $3835\:m$ \\
+ Rhone & $2200\:m$ & $3480\:m$ & $2918\:m$ & $14.4\:km^2$ & $9751\:m$ & $3256\:m$ \\
+ \end{tabular}
+ \caption{Summary characteristics of the 4 validation glaciers}
+ \label{tab:attrs}
+\end{table}
+
+\subsection{Data sources}
+
+We obtained width and overall length measurements for the 4 validation glaciers from LANDSAT 5 satellite imagery and ASTER global digital elevation models (GDEMs).
+As the LANDSAT 5 imagery has a horizontal resolution of $\pm 30$ m, we prescribe a conservative $\pm 60$ m error for glacier width measurements (error for both edges of glacier boundaries).
+ASTER GDEM data have a vertical root-mean-squared error of $\pm 15$ to $\pm 25$ m, depending on several environmental conditions (surface covering, topography, surface roughness, etc.) \cite{aster_aster_2009}.
+As our model exclusively involves mountainous snow-covered regions, we utilize the more conservative $\pm 25$ m error in calculations of bed topography and ice surface elevations.
+Bed elevation validation measurements are from modeled topographies in \cite{farinotti_method_2009, farinotti_simple_2010}, constrained using multiple GPR profiles and/or borehole depths for each glacier.
+
+The Silvretta and Gries glaciers have the best-constrained mass balances with ~50 years of published data for each \cite{wgms_fluctuations_2019}.
+In order to compare the current climatic ELA of these glaciers with our modeled ELA, we determine measured ELAs from the linearly detrended, annually-measured ELA values from 1981-2010 for both glaciers, with uncertainty calculated using a 95\% margin of error.
+The Findel Glacier has similarly well-constrained mass balance measurements from a glacier stake network \cite{wgms_fluctuations_2019}, but with a much shorter record (2005-2010).
+We compare the mean ELA over this time to the modeled ELA for Findel Glacier.
+The Rhone Glacier does not have consistent year-to-year mass balance measurements.
+Instead, we take modeled steady-state ELA estimates from air temperature correlations (1971-1990) provided in \cite{zemp_distributed_2007}.
+These ELA estimates are constrained with the few years of available stake mass balances (mean $r^2$ between measured ELA and air temperature-correlated ELA is 0.89).
+No uncertainty estimates were provided for the Rhone Glacier ELA.
+For consistency, we assume Gaussian uncertainties with bounds similar to the margins of error of the mass balances for the Silvretta, Gries, and Findel glaciers ($\pm 50$ m).
+
+\subsection{Model comparisons}
+
+\begin{figure}
+ \centering
+ \includegraphics[width=\textwidth]{Figures/Hypsometry.png}
+ \caption{Bed elevation reconstructions for the four validation sites. Yellow circles denote measured bed elevation values, black lines represent the modeled bed profile, and blue shading represents model error (2 standard deviations). Note that scales are not consistent between subplots.}
+ \label{fig:hyp}
+\end{figure}
+
+\begin{figure}
+ \centering
+ \includegraphics[width=\textwidth]{Figures/Ice.png}
+ \caption{Modeled glacier ice surfaces for the four validation glaciers. Yellow circles denote measured ice elevation values, black lines represent the mean modeled bed topography (Figure 3), blue lines represent the modeled ice surface profile, and blue shading represents model error (2 standard deviations).}
+ \label{fig:ice}
+\end{figure}
+
+\begin{figure}
+ \centering
+ \includegraphics[width=\textwidth]{Figures/Width.png}
+ \caption{Glacier width modeling for the four validation sites. Compares the overall modeled areal profile (and modeled uncertainty) with discrete measured points of each glacier’s width. Yellow circles denote width measurements for points on the glacier, black lines represent the modeled width profile, and blue shading represents model error (2 standard deviations). Note that scales are not consistent between subplots.}
+ \label{fig:width}
+\end{figure}
+
+The model results, including bed topography, ice thickness, plan-profiles, and ELAs, are summarized in Figures \ref{fig:hyp}-\ref{fig:elas} for all four validation sites.
+Most of the intermediate model outputs match measured values within error.
+Points of increased disagreement likely result mainly from local variability and the inherent smoothing caused by model fit constraints and optimization.
+Exceptions to this include the overdeepened section apparent in the Gries Glacier (Figure \ref{fig:hyp}), which represents a systemic shift in bed topography not adequately captured in the model.
+Similarly, most differences in modeled and measured ice surface (Figure \ref{fig:ice}) likely result from local variations in ice thickness of a scale finer than the input data resolution (e.g. ice crevasses), but with little effect on the final ELA estimate.
+An exception to this explanation is Findel Glacier, wherein the model appears to systemically overestimate the ice thickness, and a corresponding overestimation of the ELA by a similar magnitude (see Figures \ref{fig:ice} and \ref{fig:elas}).
+Although isolating an exact reason for this overestimation is challenging, it may be related to violations of the assumed perfect plasticity of the modeled ice or to ice thinning/downwasting due to climate disequilibrium, neither of which are accounted for in this ELA model.
+Modeled glacier width results generally closely match those recorded from satellite imagery (Figure \ref{fig:width}).
+The most noticeable exception to this is the Rhone Glacier, with a few clear outliers in the accumulation area.
+These may be related to difficulties in accurately defining the glacier boundaries in the accumulation area, or else may simply represent a more complex glacier geometry that this ELA model will not fully capture.
+Regardless, these deviations do not appear to significantly affect the final ELA estimate.
+
+\begin{figure}
+ \centering
+ \includegraphics[width=\textwidth]{Figures/ELAs.png}
+ \caption{Comparison of observed ELA values, ELA estimates using the outlined ELA model, and estimates from several empirical methods (AAR, AMA, MEG, and THAR) for the four validation glaciers. Modeled and measured estimates are within the margins of error for the 4 validation glaciers, although Gries Glacier appears to have a systemic bias. The results show a mean bias in central estimates for modeled ELAs of -14.8 m relative to measured ELAs, with no consistent direction in bias. The largest difference occurs with the Gries Glacier, with a central bias in modeled ELA of -155 m relative to the observed ELA.}
+ \label{fig:elas}
+\end{figure}
+
+Modeled ELA estimates for the four validation sites and comparisons to corresponding observed ELAs are presented in Figure \ref{fig:elas}.
+The figure additionally includes empirical estimates of ELA using the AAR, AMA, MEG, and THAR methods for each validation glacier.
+We use two AAR values(0.5 and 0.67) to show how the choice of this parameter influences the estimate.
+We assign THAR a coefficient of 0.35.
+As the MEG is equivalent to THAR=0.5, these two metrics also show the effect of varying this ratio as well.
+All four validation glaciers show agreement within errors between observed and modeled ELA values, although the central estimate for Gries glacier differs by \mytilde150 meters.
+None of the empirical methods consistently estimate the ELA within observed uncertainty for all 4 glaciers, although the AAR=0.50 and MEG iterations perform the best of the empirical methods.
+THAR (with ratio=0.35) consistently underestimates the ELA, while AMA consistently overestimates the ELA for these glaciers.
+Table \ref{tab:bias} summarizes the the mean biases over the four validation sites compared to the observed mean ELA.
+The ELA model and the AAR method (ratio=0.50) perform the best overall with mean biases of -14.8 m and 21.8 m respectively.
+
+\begin{table}[]
+ \centering
+ \begin{tabular}{c|c}
+ Method & Mean bias (m) \\
+ \hline
+ AAR=0.50 & $21.8$ \\
+ AAR=0.67 & $-58.3$ \\
+ AMA & $96.8$ \\
+ MEG & $-77.1$ \\
+ THAR=0.35 & $-314$ \\
+ ELA model & $-14.8$
+ \end{tabular}
+ \caption{ELA mean biases for different methods}
+ \label{tab:bias}
+\end{table}
+
+Likely sources of error to explain discrepancies between the observations and modeled results involve more complex considerations not accounted for with the simple ELA model.
+For instance, more complex bed topographies, differences in shading/shielding by valley walls, debris cover, and accumulation through avalanching can all affect the recorded ELA in mass balance measurements, none of which are considered in the ELA model.
+It is important to note that the model is particularly sensitive to errors in bed topography, as these values influence estimates of slope, ice thickness, and width and therefore can potentially strongly affect the final ELA estimates.
+Differences in steady-state assumptions may also be an important factor in differences between modeled and measured modern ELAs.
+The ELA model assumes steady-state conditions, whereas the annual mass balance reflects emergent climate conditions.
+Glaciers typically have either an annual mass surplus or deficit in a given year, complicating comparisons of our results to measured ELA values.
+Such a limitation, however, is inherent to all morphology-based ELA models.
+For instance, all methods significantly underestimate the ELA for Gries glacier, suggesting this glacier could be strongly out of equilibrium.
+Overall, the presented results show a high degree of confidence in the model's ability to estimate glacier ELAs (within calculated uncertainties) from relatively few geomorphic inputs, supporting the use of the presented ELA model for simple valley glaciers across a wide spectrum of bed slope geometries, glacier shapes, glacier widths, and elevation extents.
+
+The incorporation of additional variables and modeling components could address some of these limitations and improve the overall effectiveness of the model.
+For instance, the methods presented in \cite{benn_exceltm_2010, pellitero_glare_2016} also include a "shape friction factor" parameter that accounts for lateral drag in topographically-constrained glacier valleys.
+This \textit{F} factor relates the frictional lateral drag to the glacier cross-sectional area and perimeter length of the ice-topography contact \cite{pellitero_glare_2016}.
+As these required inputs are also generated during intermediary steps in our ELA model, an \textit{F} factor implementation could be incorporated into the ELA model in the future, potentially improving the ELA estimates for valley glaciers.
+In light of the sensitivity of the model to bed slope, using a non-parametric interpolation method for bed topography estimates could improve the ELA estimates, but would also require more complex modifications of the model to avoid discontinuous step changes in ice thickness and other parameters.
+
+\section{Conclusions}
+
+The model described here accurately estimates ELAs from Alpine valley glaciers of varying size, topography, and areal distribution while utilizing a small set of easily-obtained measurements.
+The model provides errors based on the physical uncertainties of model inputs, a crucial factor for determining the significance and importance of results.
+We validate the model on a set of glaciers in the Alps spanning a variety of characteristics (bed topography, size, elevation extent, etc.).
+The model performs at least as well as traditional empirical methods of ELA estimation while minimizing reliance on optimized empirical coefficients, adding uncertainty estimates, and providing insight for the sensitivity of individual glaciers to model inputs.
+Based on these validations and the physics-based nature of the model, this ELA model should serve as a robust, easily applicable, self-consistent method for ELA glacier reconstructions in varied areas, including the broader European Alps, alpine regions of the Arctic, the Southern Alps in New Zealand, and similar glaciated locations.
+The model should also be readily applicable to paleoglacier reconstructions based on preserved moraine sequences, permitting rapid and consistent comparisons of glacier changes through time and across diverse regions.
+Such studies will permit enhanced insight into the mechanisms of climate change in the past, and help us to better understand present and future changes to critical glacial and water resources in a warming world.
+
+\bibliographystyle{elsarticle-num}
+\bibliography{citations}
+
+\end{document}
\ No newline at end of file
diff --git a/Docs/ELA-model_tex/methodsx.csl b/Docs/ELA-model_tex/methodsx.csl
new file mode 100644
index 0000000..12a5f1b
--- /dev/null
+++ b/Docs/ELA-model_tex/methodsx.csl
@@ -0,0 +1,14 @@
+
+
diff --git a/Docs/ELA-model_tex/numcompress.sty b/Docs/ELA-model_tex/numcompress.sty
new file mode 100644
index 0000000..b187eff
--- /dev/null
+++ b/Docs/ELA-model_tex/numcompress.sty
@@ -0,0 +1,189 @@
+%%
+%% This is file 'numcompress'.
+%%
+%% Copyright (C) 2009-2012 River Valley Technologies
+%%
+%%
+%% This package may be distributed under the terms of the LaTeX Project
+%% Public License, as described in lppl.txt in the base LaTeX distribution.
+%% Either version 1.0 or, at your option, any later version.
+%%
+%% $Id: numcompress.sty 187 2012-08-18 09:36:35Z rishi $
+%% $URL: http://lenova.river-valley.com/svn/elsbst/trunk/numcompress.sty $
+%%
+\NeedsTeXFormat{LaTeX2e}
+\def\Fileversion$#1: #2 ${\gdef\fileversion{#2}}
+\def\Filedate$#1: #2-#3-#4 #5 #6 #7 ${\gdef\filedate{#2/#3/#4}}
+\Fileversion$Rev: 187 $
+\Filedate$LastChangedDate: 2012-08-18 15:06:35 +0530 (Sat, 18 Aug 2012) $
+\ProvidesPackage{numcompress}
+ [\filedate\space\fileversion\space numcompress (CVR)]
+\PackageWarningNoLine{numcompress}
+ {****************************************\MessageBreak
+ Package numcompress v,\fileversion\space loaded\MessageBreak
+ [Compress numbers (CVR)]\MessageBreak
+ ****************************************}
+\newif\ifdots \dotsfalse
+\newif\ifnumcompress \numcompresstrue
+
+\DeclareOption{dots}{\global\dotstrue}
+\DeclareOption{nodots}{\global\dotsfalse}
+\DeclareOption{compress}{\global\numcompresstrue}
+\DeclareOption{nocompress}{\global\numcompressfalse}
+
+\ProcessOptions
+
+\def\removeDot#1{\def\tmp{#1}%
+ \ifx\tmp\@empty\else\@removeDot#1\@nil\fi}
+
+\def\@removeDot#1\@nil{\edef\fchar{\expandafter\@car#1\@nil}%
+ \edef\rchar{\expandafter\@cdr#1!\@nil}%
+ \def\@xmltempa{.}\def\@xmltempb{!}%
+ \ifx\fchar\@xmltempb\relax\else%
+ \ifx\fchar\@xmltempa\relax\else%
+ \fchar\ignorespaces\fi\removeDot{\rchar}\fi}
+
+\def\First[#1]{\csname First#1\endcsname}
+\def\Second[#1]{\csname Second#1\endcsname}
+
+\def\parseFirstPage#1{\@tempcnta=0
+ \@tfor\@digits:=#1\do{%
+ {\global\advance\@tempcnta by 1
+ \expandafter\xdef\csname
+ First\the\@tempcnta\endcsname{\@digits}%
+ \xdef\flength{\the\@tempcnta}}}}
+
+\def\parseSecondPage#1{\@tempcnta=0
+ \@tfor\@digits:=#1\do{%
+ {\global\advance\@tempcnta by 1
+ \expandafter\xdef\csname
+ Second\the\@tempcnta\endcsname{\@digits}%
+ \xdef\llength{\the\@tempcnta}}}}
+
+\newif\ifdissimilar\dissimilarfalse
+\def\checkequal#1#2{\edef\Farg{#1}\edef\Sarg{#2}%
+ \edef\One{A}%
+ \ifcat\One\Farg \relax\else%
+ \ifdissimilar\Sarg\else%
+ \ifnum\Farg=\Sarg\relax\else\Sarg\dissimilartrue\fi\fi\fi}
+%
+\let\@@fpage\@empty
+\let\@@lpage\@empty
+\def\fpage@compress#1{%
+ \gdef\@@fpage{#1}%
+ \edef\llength{0}%
+ \parseFirstPage{#1}%
+ \ifnum\flength=\llength%
+ \gdef\@fpage{\@@fpage}%
+ \gdef\@lpage{%
+ \@ifundefined{Second1}{}{\checkequal{\First[1]}{\Second[1]}}%
+ \@ifundefined{Second2}{}{\checkequal{\First[2]}{\Second[2]}}%
+ \@ifundefined{Second3}{}{\checkequal{\First[3]}{\Second[3]}}%
+ \@ifundefined{Second4}{}{\checkequal{\First[4]}{\Second[4]}}%
+ \@ifundefined{Second5}{}{\checkequal{\First[5]}{\Second[5]}}%
+ }%
+ \else%
+ \gdef\@fpage{\@@fpage}%
+ \gdef\@lpage{\@@lpage}%
+ \fi}
+
+\def\lpage@compress#1{%
+ \gdef\@@lpage{#1}%
+ \parseSecondPage{#1}%
+ \ifnum\flength=\llength%
+ \gdef\@fpage{\@@fpage}%
+ \gdef\@lpage{%
+ \edef\One{A}%
+ \edef\xFirst{\First[1]}%
+ \edef\xSecond{\Second[1]}%
+ \ifcat\One\xSecond\relax%
+ \ifx\xFirst\xSecond%
+ \@ifundefined{Second1}{}{\checkequal{\First[1]}{\Second[1]}}%
+ \@ifundefined{Second2}{}{\checkequal{\First[2]}{\Second[2]}}%
+ \@ifundefined{Second3}{}{\checkequal{\First[3]}{\Second[3]}}%
+ \@ifundefined{Second4}{}{\checkequal{\First[4]}{\Second[4]}}%
+ \@ifundefined{Second5}{}{\checkequal{\First[5]}{\Second[5]}}%
+ \else#1\fi%
+ \else%
+ \ifx\xFirst\xSecond%
+ \@ifundefined{Second1}{}{\checkequal{\First[1]}{\Second[1]}}%
+ \@ifundefined{Second2}{}{\checkequal{\First[2]}{\Second[2]}}%
+ \@ifundefined{Second3}{}{\checkequal{\First[3]}{\Second[3]}}%
+ \@ifundefined{Second4}{}{\checkequal{\First[4]}{\Second[4]}}%
+ \@ifundefined{Second5}{}{\checkequal{\First[5]}{\Second[5]}}%
+ \else#1\fi%
+ \fi%
+ }%
+ \else
+ \gdef\@fpage{\@@fpage}%
+ \gdef\@lpage{%
+ \edef\Targ{#1}%
+ \edef\One{A}%
+ \edef\xFirst{\First[1]}%
+ \edef\xSecond{\Second[1]}%
+ \ifx\xFirst\xSecond
+ \ifcat\One\xSecond\relax\else\@@lpage\fi%
+ \else#1\fi%
+ }%
+ \fi}
+
+%\newwrite\xx
+%\immediate\openout\xx=tmpbib.tex
+\gdef\@@lpage@compress#1--#2\@nil{\lpage@compress{#1}}
+\gdef\@@@pages#1#2{\def\next{#2}%
+% \immediate\write\xx{[\the\c@NAT@ctr.]\space [#1][#2]}%
+ \fpage@compress{#1}%\ifx\next\@empty\relax\else
+ \@@lpage@compress#2\@nil%\fi
+ {\@fpage\ifx\next\@empty\relax\else
+ --\@lpage\fi}\resetall}
+
+\gdef\@@@page#1{#1\resetall}
+
+\def\mk@empty#1{\@tempcnta=1
+ \loop\ifnum\@tempcnta<6
+ \expandafter\let\csname#1\the\@tempcnta\endcsname\relax
+ \advance\@tempcnta by 1 \repeat}
+\def\resetall{\let\@lpage\@empty\let\@fpage\@empty
+ \def\flength{0}\def\llength{0}%
+ \let\@@fpage\@empty\let\@@lpage\@empty
+ \mk@empty{First}\mk@empty{Second}}
+
+
+\ifdots
+ \gdef\xfnm[#1]{\unskip\space#1}
+ \def\bibinfo#1#2{\@ifnextchar.{\@@bibinfo{#1}{#2}}{\@@@bibinfo{#1}{#2}}}
+ \def\@@@bibinfo#1#2{\def\next{#1}%
+ \def\@@@pg{pages}\def\@@@au{author}%
+ \ifx\next\@@@pg\bibpages{#2}\else
+ \ifx\next\@@@au\bibauthor{#2}\else
+ #2\fi\fi}
+ \def\@@bibinfo#1#2.{\def\next{#1}%
+ \def\@@@pg{pages}\def\@@@au{author}%
+ \ifx\next\@@@pg\bibpages{#2}.\else
+ \ifx\next\@@@au\bibauthor{#2}\else
+ #2.\fi\fi}
+\else
+ \gdef\xfnm[#1]{\unskip\space\removeDot{#1}}
+ \def\bibinfo#1#2{\def\next{#1}%
+ \def\@@@pg{pages}\def\@@@au{author}%
+ \ifx\next\@@@pg\bibpages{#2}\else
+ \ifx\next\@@@au\bibauthor{#2}\else
+ #2\fi\fi}
+\fi
+
+\ifnumcompress
+ \def\bibpages#1{\@@bibpages#1--\\\@nil}
+ \def\@@bibpages#1--#2\@nil{%
+ \ifx\\#2\relax\@@@page{#1}\else
+ \@@@pages{#1}{#2}\fi}
+ \else
+ \def\bibpages#1{#1}
+\fi
+
+\def\bibauthor#1{#1}
+
+\endinput
+
+%%
+%% End of package 'numcompress.sty'
+%%
diff --git a/Docs/Figures/mermaid-figs.md b/Docs/Figures/mermaid-figs.md
index fd87e0f..1a3cbdc 100644
--- a/Docs/Figures/mermaid-figs.md
+++ b/Docs/Figures/mermaid-figs.md
@@ -47,7 +47,7 @@ classDef default fill:#5cfcfc,stroke:#333,stroke-width:1px
"htmlLabels": false,
"useMaxWidth": false,
"rankSpacing": 15,
- "nodeSpacing": 50
+ "nodeSpacing": 100
}
}
```
diff --git a/Docs/reviewer-response.md b/Docs/reviewer-response.md
new file mode 100644
index 0000000..5e81640
--- /dev/null
+++ b/Docs/reviewer-response.md
@@ -0,0 +1,63 @@
+---
+geometry: margin=1in
+---
+
+# Response to reviewer comments
+
+## Reviewer 1
+
+> *One thing that struck me was that this is a complicated way of achieving the same results as simple techniques such as THAR and AAR. The authors note that this "ELA model is similar in simplicity to field-based geomorphic methods such as the toe-to-headwall ratio (THAR) or accumulation area ratio (AAR)". If so, the authors need to stress more convincingly why this modelling approach is superior to the simple methods. There is some justification in the text, regarding errors and uncertainties and so on and this is important but needs more prominence through the paper. One other note on the quote above is that THAR and AAR are not field-based methods. They are achievable remotely using topographic maps and satellite imagery, so you shouldn't call them field-based methods.*
+
+To address these concerns, we've added an Introduction section that better details traditional methods and highlights the key advantages of our method to address certain limitations of traditional empirical techniques (Lines 12-31).
+We've also added an additional paragraph summarizing/comparing our method to other recent ELA methods (Lines 66-81).
+The main thrust of these arguments are (1) less reliance on empirical coefficients that must be tuned using existing measured data, making their application to new regions or time periods more suspect, and (2) robust uncertainties that give greater context to the relevance of results and the sensitivity of results to specific inputs.
+Finally, we include ELA estimates from 4 various empirical methods in our validation comparisons to highlight the proposed model's abilities relative to traditional techniques (Figure 6, Table 4, Lines 222-234).
+We also specfically revised "field-based methods" to "empirical geomorphic methods" (Line 38).
+
+> *The methods of the balance equation, glacier bed modelling, ice thickness modelling and glacier width modelling are based on theoretical concepts similar to those presented in the modelling papers of Pellitero et al. 2015 and 2016. I think some discussion of the Pellitero et al models is needed as this is the most recent set of papers related to the same idea presented in your paper. I think there is plenty of room for a new modelling approach like the one proposed in your new paper, but it seems odd if no mention is made of similar modelling papers.*
+
+We have added an additional two paragraphs to the Methods section (Lines 58-81) detailing similarities and differences to other automated ELA calculations (Benn and Hulton, 2010; Pellitero et. al., 2015; and Pellitero et. al., 2016).
+
+> *I think it would also be worth noting the empirical work done on glacier ELAs. Braithwaite and Raper (2009) found that the median altitude (area-median altitude) of glaciers best described the ELAs observed on 94 modern glaciers (R2 = 0.99). How do your models compare with these empirical observations? For example, on the four modern Swiss glaciers where is the ELA positioned - is it at the area-median altitude. I would also make some recommendation as to why your approach is superior to simply using GIS to measure the mean area-median altitudes on glaciers past and present. I do like the theoretical rigour of your approach, so I am note suggesting it is unnecessary, I just want to be more convinced on why it is better to apply this more complicated approach. Or simply tell us whether it nicely complements simpler approaches such as that using a mean area-median altitude approach.*
+
+We feel the revised Intoduction section better addresses this concern, where we take greater lengths to discuss some of the previous empirical work done on ELAs (including reference to Braithwaite and Raper, 2009), as well as highlighting some of the areas our contribution can offer improvements (12-31).
+We have also included additional comparisons of the observed/modeled ELA estimates and traditional empirical techniques in our model comparison section (Figure 6, Table 4, and Lines 222-234).
+
+> *You say that the four Swiss glaciers were selected due to differences in overall shape, length and elevation extent. I've had a look at these glaciers online and they look similar sorts of glaciers to me. Can you please clarify what the differences were in the shape, length and elevation extent.*
+
+We have edited this section to better reflect that the primary rationale for glacier selection was the availability of requisite data, with the diversity in glacier characteristics an important but secondary consideration (Lines 166, 179-171).
+We have also included an additional table (Table 3) that better demonstrates the differences and similarities between the glaciers.
+
+> *Spelling error: "geoemetries" should read "geometries"*
+
+We have corrected the typo (Line 250).
+
+> *You write that "The model described here accurately estimates ELAs from a variety of glacier sizes, shapes, topographies and areal distributions". It doesn't really, since all the glaciers are similar types of valley glaciers. Perhaps it would be more precise to write "The model described here accurately estimates ELAs from Alpine valley glaciers that have a variety of sizes, shapes, topographies and areal distributions.*
+
+We have corrected the sentence to reflect that all validation was relative to Alpine valley glaciers (Line 263), while also noting in the Model Validation section that the theoretical underpinnings of the model should make it applicable to similar glaciers in other regions (Lines 172-174).
+
+## Reviewer 2
+
+> *Two references must be included within the manuscript:*
+
+>> *Benn, D.I., Hulton, N.R.J., 2010. An Excel TM spreadsheet program for reconstructing the surface profile of former mountain glaciers and ice caps. Computers & Geosciences 36, 605-610.*
+
+>> *Pellitero, R., Rea, B.R., Spagnolo, M., Bakke, J., Ivy-Ochs, S., Frew, C.R., Hughes, P., Ribolini, A., Lukas, S., Renssen, H., 2016. Glare, a GIS tool to reconstruct the 3D surface of palaeoglaciers. Computers & Geosciences 94, 77-85.*
+
+We have added an additional two paragraphs to the Methods section (Lines 58-81) detailing similarities and differences to other automated ELA calculations (Benn and Hulton, 2010; Pellitero et. al., 2015; and Pellitero et. al., 2016).
+
+> *Authors should discuss the opportunity to integrate the shape factor F in their method.*
+
+We have added an additional paragraph that discusses possible improvements to the model, including implementing an F factor and non-parametric interpolation of the bed topography (Lines 251-260).
+
+> *As a test of the ELA model's efficacy, authors compare the ELA modelled and measured for four present-day glaciers in the Swiss Alps. Readers would appreciate that the comparison be extended to other methods for ELA calculation such as AAR and AABR in order to appreciate the gap between ELA estimation from the presented physical steady-state model and empirical methods such as AAR and AABR.*
+
+The model comparisons now also include ELA estimates using AAR, THAR, AMA, and MEG in addition to the observed and modeled ELAs, highlighting some of the improvements and advantages offered by the new method (Figure 6, Table 4, Lines 222-234).
+The main results are the model does at least as well of these methods, and in some cases considerably better, while keeping the targeted advantages of the model (uncertainty estimates, sensitivity analysis, and decreased reliance on statistical coefficients).
+
+## Miscellaneous notes
+
+During revision, some of the calculated values changed slightly from the original results due to the stochastic nature of the Monte Carlo sampling.
+We have now set a specific random seed for the sampling to ensure consistent results moving forward.
+None of the changes materially impact the findings.
+The updated values are present in Figure 6, Table 4, and Lines 226-229.
diff --git a/LICENSE.md b/LICENSE.md
new file mode 100644
index 0000000..632a234
--- /dev/null
+++ b/LICENSE.md
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2020 Durban G. Keeler
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..06954d3
--- /dev/null
+++ b/README.md
@@ -0,0 +1,53 @@
+# A first-order flexible ELA model based on geomorphic constraints
+
+This page documents a flexible model to estimate steady-state equilibrium line altitudes (ELAs) of clean-ice, temperate alpine glaciers based on geomorphic inputs.
+When using this model, please cite the [original publication](https://doi.org/10.1016/j.mex.2020.101173).
+Also consider notifying me (durban.keeler@gmail.com) to help me gauge the use and interest in the project.
+
+## Background information
+
+Alpine glaciers, with their valuable combination of highly sensitive response to climate and near-global extent, are powerful tools for investigating previous and present climate changes.
+They also represent critical water resources for areas around the globe, with the potential for far-reaching effects in a warming world.
+Advancements to understand and model glacial changes and the variables influencing them are therefore relevant to many fields.
+This page contains code for a physical steady-state model for alpine glaciers relating directly to glacier mass balance (via the equilibrium line altitude) while retaining the simplicity of other morphology methods and simultaneously including error estimates.
+The model consists of custom MATLAB functions to produce ELA estimates in a user-friendly and generally-applicable method from only a limited number of glacier bed topography and glacier width measurements.
+
+## Model components
+
+The primary functions of the ELA model are located in the `src/` directory.
+The ELA model function `ELA_calc.m` requires two dataset inputs (discrete estimates of bed topography and discrete estimates of glacier width, both measured downglacier along the centerline of the glacier valley) and the number of Monte Carlo simulations to perform.
+Approximately ten quasi-equally spaced points along the length of the glacier are often sufficient, though the optimum number depends on the length and complexity of the bed topography and glacier geometry.
+To avoid issues of model extrapolation (and to automatically include the overall glacier length), both the toe and the head of the glacier should be included in these measurements.
+The ELA model input data should be provided as a MATLAB structure with four fields, as summarized in Table 1.
+Tributary glaciers, if present, should be provided as variable input arguments (formatted as MATLAB structures according to Table 1) after the number of simulations to perform.
+The `format_inputs.m` function takes .csv files of glacier bed topography and glacier width measurements and creates a properly-formatted MATLAB structure to serve as input to the ELA model.
+
+**Table 1: Required format for ELA model inputs**
+
+| Field name | Dimensions | Field description |
+|-------------|------------|-------------------------------------------------------------------------------------------|
+| `X_dist` | [N X 1] | Vector of glacier length from 0:N, where N is the total length of the glacier in meters. |
+| `Bed_pts` | [n X 2] | A matrix with positions along the glacier centerline (in meters) in the first column and corresponding bed elevation measurements (meters a.s.l.) in the second. |
+| `Ice_surf` | [n X 2] | A matrix with positions along the glacier centerline (in meters) in the first column (this should match the first column in 'Bed_pts') and corresponding ice surface elevation measurements (meters a.s.l.) in the second. |
+| `Width_pts` | [m X 2] | A matrix with positions along the glacier centerline (in meters) in the first column and corresponding glacier width measurements (meters) in the second (widths should orthogonally intersect the centerline). |
+
+In addition to the inputs, there are model parameter assumptions built into the model prescribing the assumed errors for Monte Carlo sampling.
+Updating these assumptions to better reflect specific input data is a simple matter of editing the assigned values.
+Table 2 shows a summary of these parameters and their default values.
+
+**Table 2: ELA model error assumptions**
+
+| Variable name | Default value | Variable description |
+|---------------|---------------|-----------------------------------------------------------------|
+| `zSTD` | 25 meters | Standard deviation in measured glacier bed elevation. |
+| `wSTD` | 60 meters | Standard deviation in measured glacier width values. |
+| `tau_STD` | 5E104 Pa | Standard deviation in estimated basal shear stress (used in ice thickness calculations). |
+| `rho` | 917 kg/m3 | Density of ice (used in ice thickness calculations). |
+| `g` | 9.8 m/s2 | Acceleration due to gravity (used in ice thickness calculations). |
+
+## Misc notes
+
+An example of how to use the ELA model is included in the `scripts` directory entitled `example.m`.
+This uses data from a handful of Swiss glaciers (recorded in the `Data` directory - the [original publication](https://doi.org/10.1016/j.mex.2020.101173) for full details) to demonstrate the model's use and to show some of the outputs of the model.
+The example script also utilizes a function called `shadedErrorBar` (included in the `src` directory).
+The function was developed by Rob Campbell and can be found on [his Github page](https://github.com/raacampbell/shadedErrorBar).
diff --git a/scripts/example.m b/scripts/example.m
index 95bd366..4495e7b 100644
--- a/scripts/example.m
+++ b/scripts/example.m
@@ -22,12 +22,14 @@
%% Plot results
+% Plot the modeled bed topography
figure
hold on
shadedErrorBar(vX, median(Hyp,2), std(Hyp,[],2), 'black')
scatter(glacier_data.Bed_pts(:,1), glacier_data.Bed_pts(:,2), 25, ...
'filled', 'k')
+% Plot the modeled glacier plan profile
figure
hold on
shadedErrorBar(vX, (1/2)*median(Width,2), std(Width,[],2), 'red')
@@ -37,6 +39,7 @@
scatter(glacier_data.Width_pts(:,1), -(1/2)*glacier_data.Width_pts(:,2), ...
25, 'filled', 'r')
+% Plot the modeled ice surface
figure
hold on
plot(vX, median(Hyp,2), 'black')
diff --git a/scripts/scratch.m b/scripts/scratch.m
index e4db732..0e2d868 100644
--- a/scripts/scratch.m
+++ b/scripts/scratch.m
@@ -1,14 +1,21 @@
% Script to generate ELA results and plots for validation glaciers
% Add /src directory to path
-addpath(fullfile("../src/"))
+addpath(genpath(fullfile("../src/")))
% get the folder contents
d_files = dir('../Data/');
d_list = d_files([d_files(:).isdir]==1);
data_dirs = d_list(~ismember({d_list(:).name},{'.','..'}));
-ELA_stats = zeros(length(data_dirs), 4);
+ELA_stats = table('Size', [length(data_dirs) 9], 'VariableTypes', ...
+ repmat({'double'}, 1,9), ...
+ 'VariableNames', {'AAR_50', 'AAR_67', 'AMA', 'MEG', ...
+ 'THAR_35', 'ELA_meas', 'MoE_meas', 'ELA_mod', 'MoE_mod'}, ...
+ 'RowNames', {data_dirs.name});
+
+% Set seed (for reproducibility)
+rng(0)
for i=1:length(data_dirs)
@@ -37,12 +44,24 @@
ELA_valid = detrend(elas(yr_idx)) + mean(elas(yr_idx));
ELA_MoE = 1.96*std(ELA_valid)/sqrt(length(ELA_valid));
end
-
+
+ MEG = median(mean(Hyp,2)+mean(Hx,2));
+ wid_i = mean(Width,2);
+ ice_surf = mean(Hyp,2) + mean(Hx,2);
+ AAR_idx = find(cumsum(wid_i) >= 0.50*sum(wid_i),1,'first');
+ ELA_aar50 = ice_surf(AAR_idx);
+ AAR_idx = find(cumsum(wid_i) >= 0.67*sum(wid_i),1,'first');
+ ELA_aar67 = ice_surf(AAR_idx);
+ weights = ((1/wid_i)/sum(1/wid_i))';
+ AMA = sum(weights.*ice_surf);
+ bed_topo = mean(Hyp,2);
+ thar_35 = bed_topo(end-1) + 0.35*(bed_topo(2)-bed_topo(end-1));
ELA_mu = mean(ELA_valid);
- vELA_med = median(vELA);
+ vELA_mod = median(vELA);
vELA_err = 2*std(vELA);
- stats_i = [ELA_mu ELA_MoE vELA_med vELA_err];
+ stats_i = {ELA_aar50 ELA_aar67 AMA MEG thar_35 ...
+ ELA_mu ELA_MoE vELA_mod vELA_err};
ELA_stats(i,:) = stats_i;
end
@@ -51,26 +70,70 @@
figure
hold on
-grid on
for i=1:length({data_dirs.name})
- err_meas = errorbar(2*(i-1)-0.1, ELA_stats(i,1), ELA_stats(i,2), ...
- 'blue', 'LineWidth', 3);
- pt_meas = scatter(2*(i-1)-0.1, ELA_stats(i,1), 75, 'blue', 'filled');
- err_mod = errorbar(2*(i-1)+0.1, ELA_stats(i,3), ELA_stats(i,4), ...
- 'red', 'LineWidth', 3);
- pt_mod = scatter(2*(i-1)+0.1, ELA_stats(i,3), 75, 'red', 'filled');
+ aar_67 = scatter(2*(i-1)-0.2, ELA_stats.AAR_67(i), 150, ...
+ [0.9290 0.6940 0.1250], 'filled', 'diamond');
+ aar_50 = scatter(2*(i-1)-0.1, ELA_stats.AAR_50(i), 200, ...
+ [0.9290 0.6940 0.1250], '*');
+ ama = scatter(2*(i-1), ELA_stats.AMA(i), 150, ...
+ [0.4940 0.1840 0.5560], 'filled');
+ meg = scatter(2*(i-1)+0.1, ELA_stats.MEG(i), 150, ...
+ [0.3010 0.7450 0.9330], 'filled', 'diamond');
+ thar_35 = scatter(2*(i-1)+0.2, ELA_stats.THAR_35(i), 200, ...
+ [0.3010 0.7450 0.9330], '*');
+ err_meas = errorbar(2*(i-1)-0.25, ELA_stats.ELA_meas(i), ...
+ ELA_stats.MoE_meas(i), ...
+ 'Color', [0 0.4470 0.7410], 'LineWidth', 5);
+ pt_meas = scatter(2*(i-1)-0.25, ELA_stats.ELA_meas(i), ...
+ 150, [0 0.4470 0.7410], 'filled');
+ err_mod = errorbar(2*(i-1)+0.25, ELA_stats.ELA_mod(i), ...
+ ELA_stats.MoE_mod(i), ...
+ 'Color', [0.8500 0.3250 0.0980], 'LineWidth', 5);
+ pt_mod = scatter(2*(i-1)+0.25, ELA_stats.ELA_mod(i), ...
+ 150, [0.8500 0.3250 0.0980], 'filled');
end
-legend([err_meas err_mod], {'Measured ELA', 'Modeled ELA'})
+[~,hobj] = legend([aar_67 aar_50 ama meg thar_35 err_meas err_mod], ...
+ {'AAR=0.67', 'AAR=0.50', 'AMA', 'MEG', ...
+ 'THAR=0.35', 'Observed ELA', 'Modeled ELA'}, 'FontSize', 25);
+M = findobj(hobj,'type','patch');
+set(M,'MarkerSize', sqrt(150))
ax = gca;
+set(gca, 'YGrid', 'on', 'XGrid', 'off')
+ax.GridAlpha = 0.8;
ax.XTick = 0:2:2*length(data_dirs)-1;
ax.XTickLabels = {data_dirs.name};
ax.YLabel.String = "ELA (m a.s.l.)";
ax.FontSize = 20;
hold off
+aar50_bias = mean(ELA_stats.AAR_50 - ELA_stats.ELA_meas);
+aar67_bias = mean(ELA_stats.AAR_67 - ELA_stats.ELA_meas);
+ama_bias = mean(ELA_stats.AMA - ELA_stats.ELA_meas);
+meg_bias = mean(ELA_stats.MEG - ELA_stats.ELA_meas);
+THAR_bias = mean(ELA_stats.THAR_35 - ELA_stats.ELA_meas);
+mod_bias = mean(ELA_stats.ELA_mod - ELA_stats.ELA_meas);
+% figure
+% hold on
+% plot([0 7], [0 0], 'k--')
+% scatter(1, aar50_bias, 200, [0.9290 0.6940 0.1250], '*')
+% scatter(2, aar67_bias, 150, [0.9290 0.6940 0.1250], 'filled', 'diamond')
+% scatter(3, ama_bias, 150, [0.4940 0.1840 0.5560], 'filled')
+% scatter(4, meg_bias, 150, [0.3010 0.7450 0.9330], 'filled', 'diamond')
+% scatter(5, THAR_bias, 200, [0.3010 0.7450 0.9330], '*')
+% scatter(6, mod_bias, 150, [0.8500 0.3250 0.0980], 'filled')
+% ax = gca;
+% ax.XTick = 0:7;
+% ax.XTickLabels = {'', 'AAR=0.50', 'AAR=0.67', 'AMA', ...
+% 'MEG', 'THAR=0.35', 'Modeled', ''};
+% ax.YLabel.String = "ELA mean bias (m)";
+% ax.FontSize = 20;
+% hold off
-ELA_bias = mean(ELA_stats(:,3) - ELA_stats(:,1));
+bias_tbl = table(...
+ [aar50_bias aar67_bias ama_bias meg_bias THAR_bias mod_bias]', ...
+ 'VariableNames', {'ELA mean bias (m)'}, 'RowNames', ...
+ {'AAR=0.50', 'AAR=0.67', 'AMA','MEG', 'THAR=0.35', 'Modeled'});
-gries_outERR = (ELA_stats(2,1)-ELA_stats(2,2)) - ...
- (ELA_stats(2,3) + ELA_stats(2,4));
+% gries_outERR = (ELA_stats(2,2)-ELA_stats(2,3)) - ...
+% (ELA_stats(2,4) + ELA_stats(2,5));
diff --git a/src/hyp.m b/src/hyp.m
index 6e23d82..f2f4dd3 100755
--- a/src/hyp.m
+++ b/src/hyp.m
@@ -3,4 +3,5 @@
coeff = coeffvalues(f);
Hyp = coeff(1)*exp(coeff(2)*vX) + coeff(3)*exp(coeff(4)*vX);
Hyp(Hyp<0) = 0;
+Hyp(Hyp>max(Z_pts)) = max(Z_pts);
end
\ No newline at end of file