Fracture Mechanical Analysis of Failure Processes in Antarctic Ice Shelves vom Fachbereich Maschinenbau und Verfahrenstechnik der Technischen Universität Kaiserslautern zur Verleihung des akademischen Grades Doktor-Ingenieur (Dr.-Ing.) genehmigte Dissertation von Dipl.-Ing. Carolin Plate aus Lüdenscheid Hauptreferent: Prof. Dr.-Ing. Ralf Müller Korreferenten: Prof. Dr.-Ing. Dietmar Gross Prof. Dr. Angelika Humbert Vorsitzender: Prof. Dr.-Ing. Tilmann Beck Dekan: Prof. Dr.-Ing. Christian Schindler Tag der Einreichung: 25.08.2015 Tag der mündlichen Prüfung: 01.10.2015 Kaiserslautern, 2015 D 386 Herausgeber Lehrstuhl für Technische Mechanik Technische Universität Kaiserslautern Gottlieb-Daimler-Straße Postfach 3049 67653 Kaiserslautern (cid:13)c Carolin Plate Ich danke der „Prof. Dr. Hans Georg und Liselotte Hahn Stiftung“ für die finanzielle Unterstützung bei der Drucklegung. Druck Technische Universität Kaiserslautern Hauptabteilung 5/ Bau-Technik-Energie Abteilung 5.6 Foto-Repro-Druck Alle Rechte vorbehalten, auch das des auszugsweisen Nachdrucks, der auszugswei- sen oder vollständigen Wiedergabe (Photographie, Mikroskopie), der Speicherung in Datenverarbeitungsanlagen und das der Übersetzung. ISBN 978-3-942695-11-4 Vorwort Die vorliegende Arbeit entstand während meiner Tätigkeit als wissenschaftliche Mitarbeiterin am Lehrstuhl für Technische Mechanik der Technischen Universität Kaiserslautern. Mein besonderer Dank gilt Herrn Prof. Dr.-Ing. Ralf Müller und Frau Prof. Dr. Angelika Humbert für die Ermöglichung dieser Arbeit sowie für die hervorragende Betreuung und Unterstützung sowohl von bruchmechanischer als auch von glazi- ologischer Seite. Ebenso danken möchte ich Herrn Prof. Dr.-Ing. Gross für sein Interesse an der Arbeit und für die Übernahme eines der Korreferate sowie Herrn Prof. Dr.-Ing. Beck für die Übernahme des Vorsitzes. Herzlich bedanken möchte ich mich auch bei den Kolleginnen und Kollegen für die angenehme Atmosphäre sowohl während der Arbeitszeit als auch bei den zahlreichen außeruniversitären Aktivitäten. Besonders hervorzuheben sind meine Kollegen und ehemaligen Mitbewohner Charlotte Kuhn und Markus Klassen, für die Unterstützung in allen Lebenslagen sowie meine Zimmerkollegin Julia Christmann für die intensiven wissenschaftlichen Gespräche. Zu guter Letzt möchte ich meiner Familie und meinem Freund für den Rückhalt und die fortwährende Unterstützung nicht nur während der Promotion danken. Kaiserslautern, November 2015 Carolin Plate Kurzfassung AntarktischeSchelfeisesindeinbedeutenderBestandteildesweltweitenKlimahaus- halts. Die bis zu 500m dicken schwimmenden Eisplatten sind das Bindeglied zwi- schen gegründetem Inlandeis und dem Ozean. Umgeben von Wasser auf der einen und Luft auf der anderen Seite, sind Schelfeise klimatischen Veränderungen direkt ausgesetzt. Studien der letzten Jahrzehnte zeigen, dass speziell die Schelfeise ent- langderAntarktischenHalbinselinstabilwerdenundschließlichteilweiseoderauch komplett zerbrechen. Das Zerbrechen schwimmender Schelfeisteile oder ganzer Schelfeise trägt nicht di- rekt zu einer Erhöhung des Meeresspiegels bei. Studien von Rott et al. (2002), De Angelis (2003) und Scambos (2004) zeigen jedoch, dass mit Zerbrechen des Larsen B Schelfeises, die Fließgeschwindigkeiten der dahinter liegenden Gletscher erheblich zunahm und somit verstärkt Eis aus dem Inland dem Ozean zugeführt wird. Für den sehr pessimistischen Fall eines komplett eisfreien Antarktischen Kontinents prognostizieren Fretwell et al. (2013) einen weltweiten Anstieg des Meeresspiegels um 58m. Jedoch wird auch ein Anstieg des Meeresspiegels um nur wenige Meter, wie z.B. durch das Kollabieren des Westantarktischen Eisschilds (Bamber et al., 2009b), weitreichende Konsequenzen auf die Infrastruktur entlang der Küsten haben. Es ist daher wichtig die Zusammenhänge zu verstehen, welche zu verstärktem Zerbrechen der Schelfeise führen. Im Rahmen einer bruchmechanischen Finite Elemente Analyse untersucht diese Ar- beit daher verschiedene Risssituationen in Antarktischen Schelfeisen. Basierend auf Messungen von Rist et al. (2002), werden die Materialeigenschaften des Eises für die kurze Zeitdauer des Bruchvorgangs als linear elastische angenommen. Der Last- eintrag erfolgt quasi-statisch über verschiedene Arten von Randbedingungen und Volumenkräfte. Die Bewertung der Risssituationen erfolgt über die Berechnung von Konfigurationskräften, basierend auf der grundlegenden Arbeit von Eshelby (1951). Im Zusammenspiel mit Finiten Elementen ist die Verwendung von Konfigu- rationskräften eine effektive und flexible Methode für die Berechnung von Energie- freisetzungsraten. Diese können dann in die im Ingenieurwesen für die Bewertung vonRissenüblichenSpannungsintensitätsfaktorenumgerechnetwerden. DieUnter- scheidung,obeinRisswächstodernicht,erfolgtfürModusIBelastungendurchVer- gleich mit gemessenen kritischen Werten für Spannungsintensitätsfaktoren welche zum Beispiel in Rist et al. (2002) und Christmann et al. (2015) zu finden sind. Die Lösung des zweidimensionalen linear elastischen Randwertproblems erfolgt unter Einsatz des kommerziellen Finite Elemente Programms COMSOL. Die Kon- i figurationskräfte werden in Postprocessing Routinen in MATLAB berechnet. Durch die Verwendung von Dreieckselementen mit quadratischen Ansatzfunktionen kön- nen die erwarteten Spannungssingularitäten an Rissspitzen nicht korrekt abge- bildet werden. Die gerne zur Rissbewertung herangezogene diskrete Konfigura- tionskraftanderRissspitzestelltdietatsächlichvorhandeneEnergiefreisetzungsrate daher nur ungenau dar. Eine erhebliche Verbesserung des Ergebnisses kann durch in Denzer et al. (2003) vorgeschlagenes Aufsummieren der diskreten Knotenkon- figurationskräfte in einem Bereich um die Rissspitze erzeugt werden. Bei Proble- men mit Volumenkräften, belasteten Rissflanken, räumlich verteilten Materialpa- rametern oder unsymmetrischen Geometrien und Belastungen führt die Summa- tion jedoch zu weiteren Fehlern, welche das Ergebnis verfälschen. Nur durch die Berechnung weiterer Flächenintegrale zur Berücksichtigung von Volumenkräften und räumlich verteilten Materialparametern sowie zusätzlicher Linienintegrale ent- lang der Flanken zur Einbeziehung von Druckspannungen im Inneren des Risses und gemischten Belastungsmodi, können hohe Genauigkeiten bei der Berechnung derEnergiefreisetzungsrateerzieltwerden. MitdieserImplementierungwirdderre- lative Fehler der berechneten Energiefreisetzungsrate im Vergleich zu analytischen Beispielen auf weniger als ein Promille reduziert. Die implementierten Algorithmen werden zunächst für die Analyse der Span- nungsintensitätsfaktoren vertikaler trockener Einzelrisse in Schelfeisen verwen- det. Die Datengrundlage hierfür liefern Satellitenbilder verschiedener Auf- bruchereignisse am Wilkins Schelfeis. Die Untersuchung des Einflusses der geometrischen Parameter Länge und Mächtigkeit der betrachteten idealisierten Schelfeisgeometrie sowie des Rissöffnungswinkels auf die Spannungsintensitätsfak- toren an der jeweiligen Risstiefe steht in direktem Zusammenhang mit der Be- trachtung verschiedener Randbedingungen für die Belastung der vertikalen Ränder, sowie der Abbildung des tiefenabhängigen Wasserdrucks an der Schelfeisunterseite. Hierbei zeigt sich, dass die in bisherigen semi-analytischen Studien von z.B. Van Der Veen (1998a) oder Rist et al. (2002) gerne verwendeten Spannungsrandbe- dingungen entlang der vertikalen Ränder zumindest für geringe Risstiefen ähnliche Resultate liefern, wie die in dieser Arbeit bevorzugten Verschiebungsrandbedingun- gen. Besonders hervorzuheben ist der Einfluss des Rissöffnungswinkels, welcher bei bisherigen Veröffentlichungen zu Rissen in Schelfeisen komplett vernachläs- sigt wurde. Weitere Untersuchungen mit verschiedenen Dichteprofilen, dickenab- hängigemElastizitätsmodulundverschiedenenkonstantenWertenfürdieQuerkon- traktionszahl zeigen, dass speziell die für schnelle Bruchvorgänge berechtigte An- nahmekompressiblenMaterialverhaltens(SchulsonandDuval,2009)zudeutlichen tieferen Rissen führt als das in der überwiegenden Anzahl bisheriger Studien ver- wendete inkompressible Materialverhalten. Im Folgenden wird die Analyse vertikaler Einzelrisse auf wassergefüllte Spalten ausgedehnt. Dabei werden sowohl von Schmelzwasser oder Brine gefüllte Spal- ten an der Eisoberfläche, als auch die von Meerwasser gefüllten Bodenspalten be- trachtet. Der zusätzliche Wasserdruck auf den Rissflanken sorgt bei gleicher Last ii an den äußeren Rändern für wesentlich tiefer wachsende Risse als in den bisheri- gen Rechnungen mit unbelasteten Rissflanken. Dies führt, je nach äußerer Last, und unter Annahme von kompressiblem Materialverhalten, zu einem vollständigen Durchreißen ursprünglich stabiler Risse bei einem Wasserstand von etwa der hal- ben Risstiefe. Frühere Studien mit inkompressiblem Materialverhalten von z.B. Van Der Veen (1998a) benötigen nahezu komplett gefüllte Spalten, um vollständiges Durchreißen zu erzeugen. Ein größerer, mit einem Temperatursturz einher gehender Aufbruchvorgang am Wilkins Schelfeis, motiviert die Betrachtung von Frostsprengung als möglichen Aus- löser. FürdiebruchmechanischeModellierungvonFrostsprengungwerdenzunächst mit einer thermodynamischen Simulation mögliche Eisdicken auf der Wasserober- fläche einer mit Schmelzwasser gefüllten Spalte berechnet. Der durch Phasentrans- formation im abgeschlossenen Spalteninneren entstehende Druck dient dann als Last für die bruchmechanische Auswertung. Die Simulationen zeigen, dass unter Voraussetzung einer abgeschlossenen Spalte, realistische Eisdicken zu einem voll- ständigen Durchreißen ursprünglich stabiler wassergefüllter Spalten führen. Mit einem auf Konfigurationskräfte basierenden Algorithmus zur Ermittlung der wahrscheinlichsten Rissausbreitungsrichtung untersucht die vorliegende Arbeit auch horizontal wachsende Risse in Schelfeisen. Anders als im vorherigen Fall wer- den dafür die Belastungen im kompletten Schelfeis, oder zumindest in großräumi- gen Teilen um die betrachteten Risse benötigt. Mit einem neuen Ansatz, aufbauend auf Gleichgewichtsüberlegungen für ein viskoelastisches Fluid, werden daher die aus gemessenen Geschwindigkeitsfeldern berechneten viskosen Fliessspannungen in so genannte viskose Volumenkräfte überführt. Diese dienen dann als Last für die bruchmechanische Simulation. Angewendet auf Risswachstum in der Pine Island Gletscherzunge und im Wilkins Schelfeis liefert der verwendete Algorithmus gute Übereinstimmung im Vergleich zu beobachteten Rissverläufen. Die Ergebnisse sind bemerkenswert, da keine Informationen über die räumliche Verteilung der Materi- aleigenschaften des Eises benötigt werden. Ein Nachteil der Methode ist die starke Abhängigkeit der Ergebnisse von der Qualität des gemessenen Geschwindigkeits- feldes sowie den teilweise nur unzulänglich bekannten Randbedingungen. Im Fall von hoch aufgelösten Geschwindigkeitsfeldern können aus der Berechnung der viskosen Volumenkräfte aktive Risse als Stellen mit extremen Spannungsgradienten ermittelt werden. Bei bekanntem Rissverlauf ist es darüber hinaus auch möglich, aus den Simulationsergebnissen Rückschlüsse auf mögliche Randbedingungen zu ziehen, welche dann in weiteren Analysen verwendet werden können. Für zukünftige Modellierungen von Rissen in Schelfeisen ist die Ausweitung der Algorithmen auf viskoelastisches Materialverhalten möglich. Auch eine Kopplung von Bruchmechanik und eisdynamischen Simulationen ist wünschenswert. Für eine bessere Unterscheidung von vertikal wachsenden Oberflächenspalten und horizon- tal wachsenden Rissen ist darüber hinaus eine Simulation des drei-dimensionalen Problems denkbar. iii Abstract Antarctic ice shelves are important elements of the climate system. As link between the grounded ice sheet and the ocean, the up to 500m thick floating plates are sur- rounded by water, on one side and by the atmosphere on the other side. These conditions make ice shelves very vulnerable to climate change. Recent publications indicate that, especially along the Western Antarctic Peninsula, ice shelves become unstable and eventually disintegrate. The disintegration of floating ice shelf parts or entire ice shelves does not directly contributetosealevelrise. However,studiesofRottetal.(2002),DeAngelis(2003) and Scambos (2004) show, that the disintegration of the Larsen B Ice Shelf consid- erably accelerated the flow velocities of the inflowing glaciers. This results in an increased drainage of grounded ice into the ocean. In the very pessimistic scenario of a completely ice-free Antarctic continent, Fretwell et al. (2013) predict a global sea level rise of about 58m. However, a sea level rise of only a few meters, as fore- cast by Bamber et al. (2009b) for the collapse of the West Antarctic Ice Sheet, also leads to serious consequences for the infrastructure along the continental coasts. Therefore it is of major importance to understand the processes that lead to an ac- celerated disintegration of ice shelves. ThepresentedworkanalysesdifferentfailurescenariosinAntarcticiceshelvesusing fracturemechanicalconceptstogetherwithfiniteelementsimulations. Accordingto the study of Rist et al. (2002), linear elastic material behavior during the short frac- turing process is applied. Neglecting concepts of dynamic crack propagation, the loading conditions are quasi-static using different types of boundary conditions and volume forces. The estimation of crack criticality is based on the computation of configurational forces following the approach of Eshelby (1951). In conjunction with finite ele- ments, configurationalforcesareaneffectiveandflexiblemethodfortheevaluation of energy release rates, which are then transformed into stress intensity factors, a common measure for crack criticality in engineering. The comparison of computed stress intensity factors to measured values of critical stress intensity factors allows to determine whether crack growth continues. Critical stress intensity factors for mode I opening can be found in Rist et al. (2002) or Christmann et al. (2015). The two-dimensional linear elastic boundary value problem is solved in the com- mercial finite element program COMSOL . The computation of the configurational forces follows in postprocessing routines in MATLAB . By applying triangular finite elements with standard quadratic shape functions, the expected stress singularities at the crack tip cannot be mapped. Therefore, the discrete configurational force at v thecracktipnode,frequentlyusedtocomputetheenergyreleaserate,onlyprovides an inaccurate measure for crack criticality. Denzer et al. (2003) therefore propose the summation of all nodal configurational forces in a region around the crack tip to considerably enhance the result. However, in the case of applied volume forces, spatially distributed elastic parameters, loaded crack faces or mixed crack open- ing modes, the inclusion of nodal configurational forces around the crack tip leads to further errors. Only by computing additional area integrals for applied volume forces and inhomogeneous material parameters, as well as additional line integrals alongthecrackfacestoincorporatecrackfaceloadsandmixedopeningmodes,very accurate results can be achieved. The presented algorithm shows a relative error of less than one per-mil in comparison to analytical computations. The fracture mechanical simulations first concentrate on the computation of stress intensity factors for dry single surface cracks in the ice shelf bulk. The input data forthesimulationsisbasedonsatelliteimageryofdifferentfracturescenariosinthe Wilkins Ice Shelf. The analysis of the influence of the geometric parameters length and thickness of the idealized ice shelf geometry, as well as of the crack opening an- gle is closely connected to the study of different boundary conditions for the remote loading along the vertical lateral boundaries on the one hand and for the mapping of the depth-dependent buoyancy force at the ice shelf bottom on the other hand. It appears that, at least for shallow cracks, former semi-analytical studies of e.g. Van Der Veen (1998a) or Rist et al. (2002) with stress boundary conditions along the vertical lateral boundaries are in good agreement with the results presented here, which are predominantly computed using displacement boundary conditions. The most influencing geometric parameter is the crack opening angle, which has been completely ignored in previous studies of cracks in ice shelves. Further simulations with various density profiles, depth-dependent Young’s moduli and different con- stant Poisson’s ratios show that the assumption of compressible material behavior, as motivated by e.g. Schulson and Duval (2009) for short-term fracture processes, leads to considerably deeper cracks than the use of incompressible material behav- ior, which is predominantly used in previous studies. In a next step, the examination of vertical single cracks is expanded to surface crevasses, filled by meltwater or brine, as well as to cracks at the ice shelf bot- tom. For equal lateral loading, the simulations with additional water pressure on thecrackflanksleadtosignificantlydeepercracksthanthepreviousstudieswithout crack face loads. Without additional lateral loading and under assumption of com- pressible material behavior, a water filling of about half of the crack depth results in a complete penetration of initially stable crevasses. Previous studies with incom- pressible material behavior (Van Der Veen, 1998a) and equal external loading need almost completely filled crevasses for a total penetration. A break-up event at the Wilkins Ice Shelf that coincided with a major temperature drop motivates the consideration of frost wedging as triggering mechanism for ice shelf disintegration. Prior to the fracture mechanical analysis of frost wedging, a thermodynamic simulation investigates possible ice lid thicknesses that can form on vi