An ExcelTM spreadsheet program for reconstructing the surface profile of former mountain glaciers and ice caps

ARTICLE IN PRESS Computers & Geosciences 36 (2010) 605–610

Contents lists available at ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

An ExcelTM spreadsheet program for reconstructing the surface profile of former mountain glaciers and ice caps$ Douglas I. Benn a,b,n, Nicholas R.J. Hulton c a

The University Centre in Svalbard, PO Box 156, N-9171 Longyearbyen, Norway School of Geography and Geosciences, University of St. Andrews, KY16 9AL, UK c School of Geosciences, University of Edinburgh, King’s Buildings, West Mains Road, Edinburgh EH93JW, UK b

a r t i c l e in fo

abstract

Article history: Received 11 December 2008 Received in revised form 28 August 2009 Accepted 1 September 2009

A new ExcelTM 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. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Glacier reconstruction Ice sheet surface Shape factor

1. Introduction The equilibrium-line altitudes (ELAs) of former glaciers are a widely used source of palaeoclimatic information. They provide one of the few methods for estimating palaeo-precipitation, provided independent proxies of summer temperature are available (e.g. Dahl and Nesje, 1996; Nesje and Dahl, 2003; Benn and Ballantyne, 2005). Several methods for reconstructing former glacier ELAs have been proposed (see Nesje and Dahl, 2001; Benn et al., 2005), all of which offer different advantages and disadvantages. For most situations, the most rigorous are the Accumulation Area Ratio (AAR) and Area-Altitude Balance Ratio (AABR) methods, based on well-known empirical relationships between glacier hypsometry and mass balance (Osmaston, 2005; Rea, 2009). Both methods rely on reconstructions of glacier surfaces derived from geomorphic evidence such as terminal and lateral moraines, ice-marginal meltwater channels and trimlines. Where the evidence is good, such reconstructions may be highly accurate, but where it is poorly preserved, ambiguous or lacking, former ice surfaces must be interpolated between, or extrapolated beyond, known points. This is particularly likely to be the case for former plateau ice caps, the upper surfaces of which are commonly unconstrained by ice-marginal landforms (Rea and Evans, 2005). Ice caps are especially useful as palaeoclimatic

indicators, because their ELAs are less strongly controlled by local topographic factors than are cirque and valley glaciers, and are thus more likely to provide representative data on regional precipitation patterns (Nesje and Dahl, 2001, 2003; Benn and Ballantyne, 2005). One approach to filling gaps in empirical ice surface reconstructions is to use theoretical glacier surface profiles (e.g. Schilling and Hollin, 1981; Rea and Evans, 2007). Very useful results can be obtained using simple, steady-state models which assume a ‘perfectly plastic’ ice rheology, where irrecoverable strain only occurs when the basal shear stress equals a specified yield stress (which may vary spatially). Exact and versatile solutions for the ‘perfect plasticity’ model have been available for several years, but inferior versions persist in the literature, presumably because the former have not been presented in accessible form. This paper attempts to fill this gap in the literature, by (1) outlining the main methods used for calculating equilibrium profiles of ‘perfectly plastic’ glaciers, (2) discussing some of the issues that must be considered when implementing steady-state glacier models and (3) presenting a user-friendly ExcelTM spreadsheet to allow the rapid implementation of the best method.

2. Calculating ice surface profiles $

Code available from server at http://www.iamg.org/CGEditor/index.htm Corresponding author at: The University Centre in Svalbard, PO Box 156, N-9171 Longyearbyen, Norway. Tel.: + 47 79 023367; fax: + 47 79 023301. E-mail addresses: [email protected], [email protected] (D.I. Benn). n

0098-3004/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2009.09.016

The model is built on the assumption that ice deforms in response to the driving stress (tD) (resulting from the weight and surface gradient of the ice) only when a specified yield stress (tY)

ARTICLE IN PRESS 606

D.I. Benn, N.R.J. Hulton / Computers & Geosciences 36 (2010) 605–610

is reached. If the driving stress is less than the yield stress, no ice movement can occur and the glacier will thicken and/or steepen in response to surface snow accumulation and melting, increasing tD. On the other hand, the driving stress can never exceed the yield stress, because the ice surface profile is assumed to adjust continuously to maintain the condition that tD = tY. This equality can be expressed as

tY ¼ tD ¼ rgH

@h @x

ð1Þ

where r is the density of glacier ice (  900 kg m  3), g is the gravitational acceleration (9.81 m s  2), H is the glacier thickness, h is the ice surface elevation and x is the horizontal co-ordinate, with the x-axis parallel to glacier flow, positive upglacier (Fig. 1) @h ¼ tan a @x

ð2Þ

where a is the ice surface slope. For the present purposes, it is convenient to set the origin of the x-axis at the glacier terminus, and to define a vertical z-axis with its origin at sea level. At a later stage, we will introduce a transverse horizontal (y) axis. In the literature, Eq. (1) is sometimes given with sin a in place of dh/dx, a form that is only applicable for the parallel-slab geometry and slope-parallel x-axis used in early derivations of the driving stress (see Nye, 1952; Paterson, 1994). The form of Eq. (1) given here is the correct one for a horizontal x-axis and arbitrary bed geometry (see Hooke and Le, 2005, pp. 79–81 for a full derivation). For a glacier resting on a horizontal bed, the variation of glacier thickness H along a flowline can be found by integrating Eq. (1) with respect to x. The solution is the parabolic curve Hx2 ¼ C þ

2tY x rg

ð3Þ

where C =0, following on from the condition that the ice thickness at the edge of the glacier is zero. Thus  1=2 2tY Hx ¼ x ð4Þ rg Eq. (4) can be used to obtain a quick approximation to a former glacier surface where the bed slope is small, but it is clearly inapplicable in mountainous regions where the slope of the bed can be large. For an irregular bed, an analytic solution to find the surface profile is not possible, and it is necessary to calculate it in a succession of discrete steps (numbered 1, 2, 3yi, i +1, i +2y).

We can solve Eq. (1) numerically by discretizing the surface gradient as @h h h ¼ iþ1 i Dx @x

ð5Þ

where Dx is a specified interval of distance along the x-axis. Thus, Eq. (1) can be rewritten so that the ice surface elevation at step i+1 is t  Dx Y ð6Þ hi þ 1 ¼ hi þ H i rg where H= h–B, and B is the elevation of the bed. This equation, which was introduced by Schilling and Hollin (1981) provides only an approximate solution to the problem. Since the values of tY and H used in Eq. (6) are those at step i, and may be unrepresentative of the interval between step i and i+ 1, there may be a tendency for calculated points to either overshoot or undershoot the desired value. In consequence, there is a risk that errors will accumulate, and that elements of the reconstructed ice surface may be artefacts of the numerical method. Furthermore, evaluating the ice surface elevation at step i+1 depends on knowing the ice thickness H at some initial step i. This first step cannot be at the terminus, where H=0 (division by zero), so the ice thickness at some finite distance Dx from the terminus must be specified as the starting point of the sequence. To circumvent this problem, Schilling and Hollin (1981) provided a lookup table, specifying H for a range of bed slopes. Alternatively, it may be possible to estimate H from geomorphological evidence. An exact solution to the problem was presented by Van der Veen (1999, pp. 148–150). To be fully correct, the basal shear stress and ice thickness on the right-hand side of Eq. (6) should be evaluated at step i+ 1/2, that is, they should be values calculated for the mid-point of the interval i to i +1, rather than one end of the interval. This approach yields h2i þ 1 hi þ 1 ðBi þ Bi þ 1 Þ þhi ðBi þ 1 Hi Þ

2Dx tY ¼0 rg

ð7Þ

where the overbar indicates that the yield stress is the average for the interval (see Appendix A for derivation). Van der Veen did not show how hi + 1 is to be extracted from this expression, providing an interesting exercise for numerate students of glaciology but limiting its usefulness for those without the necessary mathematical background. Although Eq. (7) looks daunting, it is a standard quadratic equation, the solution of which is straightforward. Quadratic equations have the form ay2 þ by þc ¼ 0

ð8Þ

In case of Eq. (11), the unknown y is hi + 1, and a ¼ 1; b ¼ ðBi þ Bi þ 1 Þ and 2Dx tY c ¼ hi ðBi þ 1 Hi Þ rg

ð9Þ

All the values required to calculate the coefficients a, b and c are either specified model inputs (Dx, r, g, tY, B) or are outputs from the previous (ith) step of the model (hi, Hi). The general formula for solving quadratic equations is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b 7 b2 4ac y¼ ð10Þ 2a

Fig. 1. Sketch showing definition of model geometry.

Thus, the ice surface elevation at each step can be found by inserting Eqs. (9) into Eq. (10), and iterating upglacier along the ice flowline. There is now no requirement that Hi must be nonzero, so the process can begin at the glacier terminus. Note that Eq. (10) yields the two solutions (‘roots’) of the quadratic equation. In the present case, only one of these solutions has

ARTICLE IN PRESS D.I. Benn, N.R.J. Hulton / Computers & Geosciences 36 (2010) 605–610

physical meaning (i.e. where h4B), and the other solution can be ignored. The assumption that the shear stress at the bed equals the driving stress (Eq. (1)) might be a good approximation for nonfast-flowing parts of ice sheets, but is unrealistic for ice streams, outlet glaciers, valley glaciers and other topographically controlled ice masses. In these cases, significant resistance to flow can also be provided by side-drag and differential ‘pushes and pulls’ in the direction of flow (longitudinal stress gradients). The latter are not considered in the ‘perfect plasticity’ model, but the effects of side-drag can be readily incorporated using a ‘shape factor’ (Nye, 1952). The shape factor f, represents the proportion of the driving stress that is supported by the bed, such that

tB ¼ tY ¼ f tD

ð11Þ

where tB is the basal shear stress. The remainder of the driving stress is assumed to be supported by side drag, which is a function of the shape of the confining valley or channel. Making the simplifying assumption that tB = tD at the centreline, f can be defined as

607

replacing tY with tY/f in Eqs. (4) and (9c). In case of a glacier with no side drag (e.g. an ice cap without confining walls), f=1.

3. The Profiler spreadsheet program The model described above can be implemented using the ExcelTM spreadsheet program Profiler. Two versions of the model are included. In Profiler v. 1, a single value of the yield stress is used for the whole model domain, the step length Dx is held constant, and the shape factors are not included. In Profiler v. 2, different values for the yield stress, x co-ordinate and shape factor must be entered for each step. Profiler v. 1 is useful as a teaching tool, and as for a quick ‘first look’ at a problem, while Profiler v. 2 is intended for use as a research tool. Required inputs are

ð12Þ

(1) the long profile of the bed, in the form of horizontal (x) and vertical (z) co-ordinates of a series of points along a flowline, beginning at the glacier terminus; (2) ‘target elevations’; (3) shape factors; (4) yield stress for each step.

where A is the area of the glacier cross section, H is the ice thickness and p is the ‘glacierized perimeter’ of the cross-section (Paterson, 1994). Shape factors are incorporated into the model by

Generally, it is convenient to locate each step where the long profile crosses a contour on a topographic map, as this provides an exact value for the bed elevation.

f¼

A Hp

Fig. 2. Graphical output from Profiler v. 1 and v. 2. (a) Calculated profile of glacier nourished on plateau using v. 1, with constant yield stress. (b) Calculated profile of glacier with same basal topography as (a), using v. 2 and varying yield stress. Basal stress is assumed to be zero in centre of plateau. (c) Calculated profile of glacier with steep backwall. Note unrealistic thin, steep ice cover on backwall.

ARTICLE IN PRESS 608

D.I. Benn, N.R.J. Hulton / Computers & Geosciences 36 (2010) 605–610

600 500

Elevation

400 300 200 100 0 0

200

400

600 800 1000 Distance across valley

1200

1400

1600

Fig. 3. Graphical representation of cross profile measurements required for shape factor calculations. Horizontal (y) and vertical (z) co-ordinates of former glacier bed at each point across section can be derived from topographic maps.

The ‘target elevations’ are included as an empirical guide to find the best-fit value of the yield stress. Target elevations may be based on lateral moraines, trimlines and so on. The values of yield stress fitted to the target elevations can then be used to interpolate and extrapolate the rest of the profile. Care should be taken not to accept the results at face value, as the actual basal shear stresses may have differed substantially from those reconstructed for the constrained parts of the profile. Thermal regime and bed type can vary along-flow, and these may have exerted a strong influence on the glacier profile. An additional problem exists for plateau ice caps. At the flow divide, where the ice surface is horizontal, basal shear stress is zero. If a finite yield stress is assumed, the profile will continue to rise in an unrealistic way (Fig. 2a). The basal shear stress is more likely to diminish towards zero as the divide is approached (Fig. 2b). Unfortunately, there is no simple formula for deciding the right values (and even the non-simple ways are unlikely to be accurate). This problem, however, is unlikely to be serious for small ice caps, or where ice surface elevations can be determined from trimlines on nunataks. Perhaps the best approach is to calculate alternative profiles for the summits of plateau ice caps using different stress distributions, in the hope that these bracket the ‘true’ profile. Such alternative ice surfaces have the additional advantage of allowing errors associated with the reconstruction process to be estimated. Another issue with implementation of the program arises where the former glacier had a steep backwall. Since the model is based on the assumption of a constant basal stress, calculated glacier thicknesses will always be finite. This means the glacier surface cannot intercept the backwall, as real glaciers can. Instead, the program produces a very thin, steep glacier on the backwall (Fig. 2c). It is possible to remedy this problem numerically, but for simplicity this is not included in the present version of Profiler. The upper limit of the former glacier will be very closely approximated by the sharp inflection in the reconstructed surface, and can be readily identified visually. A third sheet, shape factor calculator will calculate shape factors from cross-profile distance and bed elevation data (y and z coordinates), which can be easily derived from topographic maps (Fig. 3). The elevation of the end points of the cross profile (which should be equal) are only known with certainty where clear geomorphic evidence exists. Where such evidence is lacking, shape factors must be calculated using estimated glacier surface elevations. If these are found to differ too greatly from the

calculated long profile, shape factors can be recalculated using calculated elevations as input, but further iterations of this procedure are unlikely to be necessary.

4. Conclusions Profiler provides a robust and flexible tool for reconstructing former glacier long profiles, using an exact solution of a ‘perfectly plastic’ glacier model. Profiler v. 1 can be used as a teaching tool, to introduce students to the principles of glacier reconstruction, and allow them to experiment with varying yield stress. Profiler v. 2 is designed for use as a research tool, and for teaching more advanced students. The availability of a user-friendly program for calculating glacier long profiles will, it is hoped, facilitate the reconstruction of palaeo-glaciers by geomorphologists and geochronologists.

Acknowledgements Careful reviews by Simon Carr and Brice Rea added substantially to the clarity and consistency of the paper. We are also grateful to Felix Ng for discussions about certain aspects of model behaviour.

Appendix A Derivation of Eq. (7) Starting with Eq. (6) t  Dx y hi þ 1 ¼ hi þ H i rg

ðA1Þ

we need to evaluate H at the position i+ 1/2 for the approximation to the gradient term – dh/dx in Eq. (2) – to be properly defined for the interval from hi to hi + 1. Thus we can put Hi þ 1=2 ¼

ðHi þ Hi þ 1 Þ 2

Substituting back into (A1) we now get   ty Dx hi þ 1 ¼ hi þ ðHi þHi þ 1 Þ=2 i rg

ðA2Þ

ðA3Þ

ARTICLE IN PRESS D.I. Benn, N.R.J. Hulton / Computers & Geosciences 36 (2010) 605–610

where the overbar on ty denotes that the yield/driving stress is averaged over the same interval hi to hi + 1. This then gives us ðhi þ 1 hi ÞðHi þ Hi þ 1 Þ ¼

2ty Dx rg

ðA4Þ

making the further substitution of Hi + 1 = hi + 1 Bi + 1, we get ðhi þ 1 hi ÞðHi þ hi þ 1 Bi þ 1 Þ ¼

2ty Dx ri g

ðA5Þ

expanding, we then get 2

hi þ 1 Hi þ hi þ 1 hi þ 1 Bi þ 1 hi Hi hi hi þ 1 þ hi Bi þ 1 ¼

2ty Dx ri g

ðA6Þ

Lastly we further substitute the first term as hi + 1Hi = hi + 1(hi  Bi) to give 2

hi hi þ 1 hi þ 1 Bi þ hi þ 1 hi þ 1 Bi þ 1 hi Hi hi hi þ 1 þ hi Bi þ 1 ¼

2ty Dx ri g ðA7Þ

with a simple rearrangement of terms yields Eq. (7): h2i þ 1 hi þ 1 ðBi þ Bi þ 1 Þ þ hi ðBi þ 1 Hi Þ

2 Dxty ¼0 ri g

ðA8Þ

Appendix B B.1. Profiler ExcelTM spreadsheet: additional explanation These notes expand on the operation of the Profiler Excel Spreadsheet described above.

B.2. Worksheet ‘Profiler v. 1’ The glacier is represented by calculating new surface elevations at discrete x locations along the horizontal flow line of the glacier/ ice sheet. Each location in x is represented by a row of the data in the main data block in rows 9 and below. Thus the ice sheet profile evolves from the margin to the interior ‘down the page’, where a given row is in fact equivalent to the ‘i+ 1th’ location when evolving the surface in the horizontal as given in Eq. (10) of the paper. Thus the preceding row (above), refers to the ‘ith’ location in equations that require evolution forwards. A few constants are held in a separate data header. Columns B:F hold data describing the system, Columns Q and R are used additionally as calculation space. Columns G:P are unused except to accommodate the graph. Values that need to be entered by the user are designated by bold red text, those that are calculated by black text.

B.3. Explanation of data B.3.1. Data header (Rows 2:5, Column C). Row: 2: Ice density: The density of ice (ri) in kg m  3. 3: g: Gravitational acceleration (g) in ms  2. 4: Yield stress: Yield stress of ice (ty) in Pa (Pascals), as given in the assumptions of perfect plasticity. This is the principal variable to change to provide different surface profiles. In this version, yield stress is the same across the profile (c.f. Profiler v. 2 below). 5: Step length: The distance in metres between discrete calculation points in x. In this version, the step length is constant along the profile.

609

B.4. Main data block (Columns B:F, Q:R) B: Distance from terminus (m): Calculated incrementally using the starting distance (row 9), and adding in the step length (cell C5) in each horizontal step. C: Bed elevation (m): This is data that must be entered by the user to provide a bed elevation at each point along the profile. D: Target elevations (m): These are also entered by the user and are used for information/display only. They are not used directly in any calculations. They appear on the graph as a heavy blue line. The idea is that the model profile can be adjusted towards the line showing these values by modifying the global yield stress in cell C4. E: Ice surface (m): In Row 9, the ice surface is set to the bed elevation ( = C9), thus assuming zero ice thickness. In Rows 10 and below the elevation is progressively evolved based on the ice surface and thickness in the preceding row (above, i.e. toward the margin), and the bed elevation for this location/row. These calculations are done here and in columns, F, Q and R which evaluate different parts of the overall calculation required as successive dependencies. In this column, Eq. (10) in the paper is evaluated, thus, e.g.: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b 7 b2 4ac ð10Þ y¼ 2a where y is equivalent to surface elevation, hi + 1, evaluates for Row 10 to the ExcelTM equation ðQ 10 þ ððQ 102 Þð4 R10ÞÞ0:5 Þ 2 The value b in Eq. (9), is calculated in Column Q, and c in Eq. (10), from Column R. In Eq. (9), a =1 so does not need to be evaluated. As explained in the paper, Eq. (10) solves for two roots, but the evaluation here is for the one that is physically meaningful (i.e. the ice surface is above the bed). F: Ice thickness (m). This is calculated simply as the difference between the ice surface (Column E) and the bed (Column C). Thus e.g., cell F10 evaluates to E10–C10. Q: This evaluates the quotient ‘b’ using Eq. (9) in the paper b=  (B + Bi + 1) where B is the bed elevation. Thus we add the bed elevation for ‘this’ location to the one at the next (i+ 1th) location, which gives the equivalent ExcelTM equation as =  (C9+C10) R: This column evaluates the quotient ‘c’ using Eq. (9) in the paper 2Dx ty C ¼ hi ðBi þ 1 Hi Þ ri g and is in fact the most complex part of the calculation. This evaluates to, e.g., row 10: =(E9(C10–F9))  ((2(B10–B9)$C$4)/($C$2$C$3)) where the values held in $C$2, $C$3 and $C$4 are respectively the constants, r(ice density), g (gravity) and ty described in the data header above. Dx (step length), is calculated explicitly as the difference between horizontal locations held in Column B (here, B10–B9). Since ‘this’ row (here row 10) denotes the position ‘i+1’ then, the position ‘i’ is designated by row 9. Thus hi is given here by cell reference E9, Bi + 1 by cell C10 and Hi by cell F9. B.5. Worksheet ‘Profiler v. 2’: The overall layout of this worksheet is more or less identical to that of Profiler v. 1. The principal differences are

 Yield stress is defined separately for each horizontal location (in Column G). Thus Row 4 in the data header is empty and the

ARTICLE IN PRESS 610



D.I. Benn, N.R.J. Hulton / Computers & Geosciences 36 (2010) 605–610

user must enter a yield stress value (in Pa) for each row in Column G. A user-defined shape factor (f) is typed in for each location (Column H). This modifies the yield stress as given in Eq. (14) in the paper and explained in the text.

As a consequence, the equation evaluating the quotient c in Column Q is modified accordingly.

Appendix C. Supplementary material Supplementary material associated with this article can be found in the online version at doi:10.1016/j.cageo.2009.09.016.

References Benn, D.I., Ballantyne, C.K., 2005. Palaeoclimatic reconstruction from Loch Lomond Readvance glaciers in the West Drumochter Hills Scotland. Journal of Quaternary Science 20, 577–592. Benn, D.I., Owen, L.A., Osmaston, H., Seltzer, G.O., Porter, S.C., Mark, B., 2005. Reconstruction of equilibrium line altitudes for tropical and sub-tropical glaciers. Quaternary International 138-139, 8–21.

Dahl, S.O., Nesje, A., 1996. A new approach to calculating Holocene winter precipitation by combining glacier equilibrium line altitudes and pine-tree limits: a case study from Hardangerjøkulen, central southern Norway. The Holocene 6, 381–398. Hooke, R., Le, B., 2005. Principles of Glacier Mechanics 2nd ed. Cambridge University Press, Cambridge 429 pp. Nesje, A., Dahl, S.O., 2001. Glaciers and Environmental Change.. Hodder Arnold, London 216 pp. Nesje, A., Dahl, S.O., 2003. Glaciers as indicators of Holocene climate change. In: Mackay, A., Battarbee, R., Birks, J., Oldfield, F. (Eds.), Global Change in the Holocene. Hodder Arnold, London, pp. 264–280. Nye, J.F., 1952. The mechanics of glacier flow. Journal of Glaciology 2, 82–93. Osmaston, H., 2005. Estimates of glacier equilibrium-line altitudes by the area  altitude, area  altitude balance ratio and the area  altitude balance index methods and their validation. Quaternary International 138-139, 22–31. Paterson, W.S.B., 1994. The Physics of Glaciers 3rd ed. Butterworth-Heinemann, London 480 pp. Rea, B.R., 2009. Defining modern day area-altitude balance ratios (AABRs) and their use in glacier-climate reconstructions. Quaternary Science Reviews 28, 237–248. Rea, B.R., Evans, D.J.A., 2005. Plateau icecap landsystems. In: Evans, D.J.A. (Ed.), Glacial Landsystems.. Hodder Arnold, London, pp. 407–431. Rea, B., Evans, D.J.A., 2007. Quantifying climate and glacier mass balance in north Norway during the Younger Dryas. Palaeogeography, Palaeoclimatology, Palaeoecology 246, 307–330. Schilling, D.H., Hollin, J.T., 1981. Numerical reconstructions of valley glaciers and small ice caps. In: Denton, G.H., Hughes, T.J. (Eds.), The Last Great Ice Sheets. Wiley, New York, pp. 207–220. Van der Veen, C.J., 1999. Fundamentals of Glacier Dynamics. Balkema, Rotterdam 462 pp.