Repositorio Institucional de Documentos

2014 33 David Chinarro Vadillo System Engineering Applied to Fuenmayor Karst Aquifer (San Julián de Banzo, Huesca) and Collins Glacier (King George Island, Antarctica) Departamento Director/es Informática e Ingeniería de Sistemas Villarroel Salcedo, José Luis Cuchi Oterino, José Antonio Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es David Chinarro Vadillo SYSTEM ENGINEERING APPLIED TO FUENMAYOR KARST AQUIFER (SAN JULIÁN DE BANZO, HUESCA) AND COLLINS GLACIER (KING GEORGE ISLAND, ANTARCTICA) Director/es Informática e Ingeniería de Sistemas Villarroel Salcedo, José Luis Cuchi Oterino, José Antonio Tesis Doctoral Autor 2014 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA System Engineering Applied to Fuenmayor Karst Aquifer (San Julián de Banzo, Huesca) and Collins Glacier (King George Island, Antarctica) David Chinarro Vadillo PhD. Dissertation Doctoral advisors: Dr. José Luis Villarroel Salcedo Dr. José Antonio Cuchí Oterino Dissertation submitted for the PhD. degree of Computing and Systems Engineering University of Zaragoza Zaragoza, October 2013 ... to my wife, for her devoted support and infinite love ... to my children, to encourage them in the adventure of acquiring and sharing new skills ... to my mother, who instilled me the importance of hard study from early age Resumen LAingeniería de sistemas, definida generalmente como arte y ciencia de crear soluciones integrales a problemas complejos, se aplica en el presente documento a dos sistemas naturales, a saber, un sistema acuífero kárstico y un sistema glaciar, desde una perspectiva hidrológica. Las técnicas de identificación, desarrolladas típicamente en ingeniería para representar sistemas artificiales por medio de modelos lineales y no lineales, pueden aplicarse en el estudio de los sistemas naturales donde se producen fenómenos de acoplamiento entre el clima y la hidrosfera. Los métodos evolucionan para afrontar nuevos campos de identificación donde se requieren estrategias para encontrar el modelo idóneo adaptado a las peculiaridades del sistema. En este sentido, se han considerado especialmente las herramientas basadas en la transformada wavelet utilizadas en la preparación de series temporales, suavizado de señales, análisis espectral, correlación cruzada y predicción, entre otros. Bajo este enfoque, una aplicación a mencionar entre las tratadas en esta tesis, es la determinación analítica del núcleo efectivo estacional (SEC) a través del estudio de la coherencia wavelet entre temperatura del aire y la descarga del glaciar, que establece un conjunto de períodos de muestreo aceptablemente coherentes, a partir del cual se crearán los modelos del sistema glacial. El estudio está dirigido específicamente a estimar la influencia de la precipitación sobre la descarga del acuífero kárstico de Fuenmayor, en San Julián de Banzo, Huesca, España. De la misma manera, se ocupa de las consecuencias de la temperatura del aire en la fusión del hielo glaciar, que se manifiesta en la corriente de drenaje del glaciar Collins, isla King George, Antártida. En el proceso de identificación paramétrica y no paramétrica se buscan los modelos que mejor representen la dinámica interna del sistema. Eso conduce a pruebas iterativas, donde se van creado modelos que se verifican sistemáticamente con los datos reales del muestreo, de acuerdo a un criterio de eficiencia dado. La solución mejor valorada según los resultados obtenidos en los casos tratados apuntan a estructuras de modelos en bloques. Esta tesis significa una exposición formal de la metodología de identificación de sistemas propios de la ingeniería en el contexto de los sistemas naturales, que mejoran los resultados obtenidos en muchos casos de la hidrología kárstica que comúnmente usaban métodos ad hoc ocasionales de carácter estadístico; así mismo, los enfoques propuestos en los casos de glaciología con el análisis wavelet y los modelos orientados a datos raramente considerados en la literatura, revelan información esencial ante la imposibilidad de precisar la totalidad de la física que rige el sistema. Notables resultados se derivan en la caracterización de la respuesta del manantial de Fuenmayor y su correlación con la precipitación, desde la perspectiva de un sistema lineal, que se complementa con los métodos de identificación basados en técnicas no lineales. Así mismo, la implementación del modelo para el glaciar Collins, obtenido también mediante métodos de identificación de caja negra, puede revelar una inestabilidad de los límites de los periodos activos de la descarga, y consecuentemente la variabilidad en la tendencia actual en el cambio climático global. Palabras clave. Hidrología, karst, glaciología, glaciar, wavelet, coherencia, identificación de sistemas, ingeniería de sistemas, sistema lineal, sistema no lineal 3 Contents Contents 4.3.3 Infrastructure and facilities ................................................116 4.4 Analysis of signals .............................................................118 4.4.1 Sampling period ..........................................................118 4.4.2 Time series preparation ...................................................119 4.4.3 Linear analysis. Simple correlation and spectral analysis ....................120 4.4.4 Wavelet analysis .........................................................121 4.4.4.1 Wavelet power spectrum of rainfall signal ...........................121 4.4.4.2 Wavelet power spectrum of discharge signal ........................123 4.5 System identification ...........................................................125 4.5.1 Nonparametric system identification .......................................125 4.5.1.1 Cross-correlation. Impulse response ...............................126 4.5.1.2 Wavelet Coherence ...............................................127 4.5.2 Parametric system identification ...........................................129 4.5.2.1 Identification of the impulse response ..............................129 4.5.2.2 Identification of transfer functions ..................................129 4.6 Nonlinear Identification .........................................................132 4.7 Conclusion ....................................................................133 5 Collins glacier 135 5.1 Introduction ....................................................................135 5.2 Geographical framework ........................................................137 5.2.1 King George Island .......................................................137 5.2.2 KGI Climate..............................................................139 5.3 Collins glacier..................................................................140 5.3.1 Hydrological overview ....................................................140 5.3.2 Instruments and facilities ..................................................142 5.3.2.1 Discharge volume measurement ...................................142 5.3.2.2 Weather measurements ...........................................142 5.4 Analysis of time series ..........................................................143 5.4.1 Time series preparation ...................................................144 5.4.2 Discharge ...............................................................144 5.4.3 Air temperature ..........................................................145 5.4.4 The involvement of air temperature in the discharge .........................145 5.4.5 Time series power spectrum ..............................................146 5.4.6 Periods of the annual cycle................................................148 5.4.7 Seasonal Effective Core calculus ..........................................150 5.5 System identification ...........................................................153 5.5.1 Linear model by cycles ...................................................154 5.5.2 Nonlinear model by cycles ................................................155 5.5.3 Glacier Coalescent Model (GCM) ..........................................157 5.6 Conclusion ....................................................................158 6 Final conclusions 161 6.1 Introduction ....................................................................161 6.2 Synthesis......................................................................162 6.3 Summary of the main chapters ..................................................163 6.3.1 Analysis and Identification of Fuenmayor aquifer ............................164 6.3.2 Analysis and Identification of Collins glacier ................................165 6.4 Findings .......................................................................165 6.5 Contributions ..................................................................166 6.6 Discussions about results .......................................................167 6.7 Perspectives ...................................................................169 6.7.1 Suggestion for onward research lines.......................................169 6.7.2 Final Remarks ...........................................................170 References .........................................................................171 Appendices 187 11 Contents Contents A Abstracts 187 A.1 International Journals and books ................................................187 A.2 International Congresses .......................................................189 A.3 Other works related to this dissertation...........................................190 B Graphics 195 C Glossary 201 12 List of Figures 1.1 Estimates referred to volumes of glaciers and ice caps (left), and groundwater (right). Both are the major topics to target mathematical models of this thesis. Source: Rekacewicz (2009) .................................................... 24 1.2 Checkland’s taxonomy (Zexian and Xuhui, 2010): natural system, designed physical system (engineered system), designed abstract system and human activity system. Climate and karst (or glacier) are coupled systems in the context of system theory. Natural systems can leverage the tools from systems engineering 28 2.1 Morlet’s «wavelet daughters» ψs,τ(t), that are dilated and translated functions derived from the «mother wavelet» ψ(t), according with Eq.2.8 ................... 44 2.2 Big-O complexity for Fourier Transform, FFT and DWT ........................... 51 2.3 Hammerstein-Wiener model. ................................................... 56 2.4 Identification problem .......................................................... 57 2.5 Stages and components in the system identification problem ..................... 58 2.6 System model structures. ...................................................... 67 3.1 The essential components of the karst aquifer (White, 2003) ...................... 78 3.2 Interface of a glacier system. The upper boxes indicate the state variables. Lowwer boxes indicate the processes involved in interactions oriented according the arrows. 90 3.3 Earth’s annual and global mean energy balance from Kiehl and Trenberth (1997). . . 91 3.4 The hydrology development and melting process in a glacier, modified from Price et al. (1979)................................................................... 94 3.5 Impact weather variables on glacier discharge, categorized by melting models and runoff models. ................................................................ 95 3.6 Diagram illustrating that the exergy component due to transfer accompanying heat is the responsible of governing the melting process and consequently the discharge 99 4.1 Huesca province containing Parque Natural de la Sierra y Cañones de Guara. The green triangle is the location of Fuenmayor spring. Map from Cuchí and Setrini (2004). .......................................................................109 4.2 An overview of San Julián de Banzo town and basin of Fuenmayor spring. ........111 4.3 Aerial view of the source area of Fuenmayor spring, marked by a red circle. In both side of the stream, in the left corner of the picture, are the districts of San Julian de Banzo (Suso to the left, and Yuso to the right). Courtesy of Google-Map. .......112 4.4 Geological scheme of the study area, according to Millán (1996) ..................113 4.5 Schema of two structural hydrogeological units: Fuenmayor (center), with its hypothetical basin limits, and Dos Caños (south) (Cuchí and Villarroel, 2002) on geological section. ............................................................114 4.6 Datalogger used in the monitoring station, at Fuenmayor spring ...................117 4.7 Signals used for the analysis of the Fuenmayor spring: (red) effective rainfall; (blue) measured discharge. ..........................................................119 4.8 Periodograms for the temporal series of Fuenmayor spring: (up) effective rainfall; (down) measured discharge. ...................................................120 13 List of Figures List of Figures 4.9 Discharge autocorrelation and spectrum for Fuenmayor spring: (a) Autocorrelation; (b) Spectral density............................................................121 4.10 Up: Continuous wavelet power spectrum of rainfall in Fuenmayor basin. Down: Continuous wavelet power spectrum of discharge time series of the Fuenmayor spring. The y axis is the wavelet scale (hours). The x axis is time (months). Curved lines on either side indicate the cone of influence where edge effects become important. The thick black line indicates significance level (established at 5%) for power spectrum. Spectral strength is shown by colors ranging from deep blue (weak) to deep red (strong). The thick black contour designates the 5% significance level against red noise. The dashed line is the cone of influence. The standardized rainfall series has an AR1 coefficient of 0.02. .......................122 4.11 Wavelet spectral power of the rainfall (up) and the discharge (down) in Fuenmayor spring zoomed to display details in high frequency bands. ........................124 4.12 cross-correlation between effective rainfall and discharge in Fuenmayor spring .....125 4.13 Kernel estimation of Fuenmayor spring based on cross-correlation: (a) without prewhitening and (b) with prewhitening. .........................................126 4.14 Linear coherence function of Fuenmayor spring..................................127 4.15 Coherence spectrum between rainfall and discharge in Fuenmayor spring. The arrows represent the phase angles between the effective rainfall and the discharge in the respective frequency band ...............................................128 4.16 Kernel estimation of Fuenmayor karst system based on: (a) Wiener-Hopf summation equation and (b) error minimization. .................................130 4.17 Impulse response and simulated runoff in the parametric identification of Fuenmayor spring . ............................................................131 4.18 Prediction of discharge in the first half of year 2010, obtained by parametric identification of the transfer function. ............................................131 4.19 Verification of model M3by Hammerstein-Wiener structure for Fuenmayor spring . . 133 5.1 Map of King George, South Shetland archipelago. Antarctic Peninsula. ...........137 5.2 Location of the catchment plot at Maxwell Bay, next to the Uruguayan base Artigas.(Source: Institut für Physische Geogrphie, Universität Freiburg, Germany. Laboratório de Pesquisas Antárticas e Glaciológicas, Universidade Federal do Rio Grande do Sul, Brazil ) ........................................................138 5.3 Topographic sketch of ice thickness at Low dome also known as Bellinghausen dome, Collins glacier. .........................................................141 5.4 Bellingshausen monthly temperature statistics: this plot shows the mean, quartiles and range for each month ......................................................143 5.5 Trend of air temperature in the period considered, with a rate of one degree per decade. ......................................................................146 5.6 Measurements of the runoff from glacier Collins at station CPE-KG-62S. a) Time series of the air temperature. b). Idem of the discharge. The third cycle data is lost for failure of sensor. ...........................................................147 5.7 Discharge wavelet spectrum at the Collins glacier. The vertical axis has a logarithmic scale in days (hours). The abscissa axis represents the time line labeled quarterly. Values of the spectral power are depicted by a colour map. The cone of influence is also represented to indicate the useful area of the spectrum without influence of the edge effects.............................................148 5.8 Air temperature wavelet spectrum at the Collins glacier. ..........................149 5.9 An approach of inactive, effective and transition periods boundaries of discharge and air temperature signals for one annual cycle. ................................149 5.10 Wavelet spectrum coherence 2005-2006 summer. The vertical colored lines indicate the limits of the head and tail of the SEC. ...............................151 5.11 Time series, wavelet coherence spectrum, and coherence levels in the 2005-2006 cycle of Collins glacier. ........................................................153 5.12 The Hammerstein block is a dead zone with a breakpoint in the temperature T0.. . . . 155 14 List of Figures List of Figures 5.13 Two different step responses. M6has a typical overdamped response (M1,M5and M8have similar efficiency). However, M2has a non negligible overshoot (M4,M7 and M10 have similar efficiency). ................................................156 B.1 3D representation of Continuous Wavelet Transform of discharge signal in Fuenmayor aquifer ............................................................195 B.2 Precipitation in the Guara area with downscaling to Fuenmayor spring station. According to the Technical Paper VI (2008) of Intergovernmental Panel on Climate Change (IPCC), the best-estimate in the global surface temperature from 1906 to 2005 is a warming of 0.74 grades C, with a more rapid warming trend over the past 50 years. The temperature increasing impacts the hydrologic cycle. The gray area in the figure represents the variability of the predicted data. ......................196 B.3 Forecasting of the discharge in Fuenmayor spring. Estimation of future trend applying Global Climate Models, downscaling methods for Aragon (Ribalaygua et al., 2013) and Hammerstein-Weiner model of this thesis (Sec. 4.6)..............196 B.4 The discharge (a) and the air temperature (b) wavelet spectrum at the glacier Collins. The vertical axis shows a logarithmic scale from 25up to 214 hours. Numbers outside the parentheses indicate the equivalent in days. The abscissa axis represents the timeline labeled quarterly, although the calculation of spectral power is held every hour. The abscissa axis represents the timeline labeled quarterly, although the calculation of spectral power is held every hour. Values of the spectral power are depicted by a distribution of color hue depending on the moment in time series and the frequency (or scale).The cone of influence indicates the useful area of the spectrum without influence of the edge effect. ...............197 B.5 Trend of the discharge of Collins glacier in the period considered (green line), after applying a wavelet filter (red line) to smooth the signal. The trend reaches 0.8 m3/s at the end of the period. .......................................................198 B.6 Wavelet spectrum coherence between air temperature and discharge of Collins preglaciar stream. ............................................................198 B.7 Collins glacier models. Verification of Model 8 (Hammerstein-Wiener) with messured data during cycle 2 (2002-2003 season) ...............................199 B.8 Collins glacier models. Verification of Model 7 (Hammerstein-Wiener) with messured data during cycle 4 (2004-2005 season) ...............................200 15 List of Tables 2.1 History of wavelets. Source: Daubechies et al. (2001) and information gathered from wavelet literature ............................................................... 41 2.1 (Continuation) .................................................................. 42 2.2 Requirements of a wavelet function ψ∈L2(Rd)(ˆ ψis the Fourier transform of the wavelet function ψ)............................................................ 43 3.1 Coefficient C in the Eagleman’s expression ....................................... 82 3.2 Classification of glaciers with the main characteristics, according to Miller (1976) ... 92 4.1 Literature review of studies about Fuenmayor spring. .............................114 4.1 (Continuation) ..................................................................115 4.1 (Continuation) ..................................................................116 5.1 Limits of the SEC for all cycles for a coherence level of 0.8. ........................152 5.2 OE linear model efficiency. M1, ... M9 refer to the models created on known data in each year ......................................................................154 5.3 Hammerstein-Wiener model efficiency. M1, ...,M10 are models based on recorded data of years 2001,..., 2011, respectively .........................................155 5.4 Average temperature in Celsius degrees. .........................................156 5.5 Characteristic parameters of the step response year by year and linear transfer function. .......................................................................157 5.6 Eficiencies of overdamped and overshoot models. ................................158 17 List of symbols Lists of symbols used in this paper with a brief description. µS/cm 25 ◦C Micro-Siemens per centimetre. Conductivity unit. θParameters vector un(t)Discrete time series for the input of system yn(t)Discrete time series for the output of system ˆ yn(t)Estimate value of yn(t) f*(x) Complex conjugate of function f(x) F{u(t)}Fourier transform of u(t) f F{u(t)}Discrete Time Fourier Transform of u(t) b f(ω)Fourier transform of the function f in the frequency domain e f(τ,ω)Wavelet transform of the function f in the time-frequency space Gψ,s,τ{u(t)}Wavelet power spectrum of u(t) with wavelet family ψ H W f1,[nl],f2,[nl],f3,[nl]Hammerstein-Wiener model from u(t) and y(t) data, and blocks f1,[nl],f2,[nl],f3,[nl] L2(R)Set of all measurable functions in the Hilbert space that are square integrable (f⊗g)(t)Convolution of two function f(t) and g(t) (f∗g)(t)Correlation of two function f(t) and g(t) L{y(t)}Laplace transforms of y(t) Wψ,s,τ{u(t)}Continuous Wavelet Transform of u(t) with wavelet family ψ Pψ,s,τ{u(t)}Wavelet power spectrum of u(t) with wavelet family ψ Cψ,s,τ{u(t),y(t)}Wavelet coherence spectrum Γψ[f(t),g(t)]τCoherence Average Function σ2 yVariance of y(n) Ck yAutocovariance of y(n)with lag k ryAutocorrelation of y(n) 19 INTRODUCTION 1.3 Basic concepts 1.3 Basic concepts System is a collection of interrelated elements that form a whole and with general properties of the whole rather than of the individual elements (Bertalanffy, 1968). The word «system» derives from the Greek «synhistanai» (σνστηµα) which means «to place together». In the ample sense, the term «system» may mean an engineered system, a natural system, a social system, or all three. Systems science provides methods to address complex problems, which enable researchers to examine the dynamic interrelationships of variables at multiple levels of analysis, and study the impact on the behavior of the system as a whole over time. Moreover, simulation modeling can be used to generate forecasting, allowing decision makers to simulate the impact of alternative solutions before carrying them out (Sterman, 1994). Checkland classified four types of systems about the real world after considering Boulding’s hierarchy and Jordan’s taxonomy. The four systems include natural system, engineered system,abstract system and human activity system (Zexian and Xuhui, 2010). Hence, the concept of system serves also to identify those manifestations of natural phenomena and processes with complex relationships among them. Anything that can belong to the Earth domain is a natural system, component of a natural system, or an aggregate formed by natural systems. In the searching the appropriate method and methodology for specific situations, according with Zexian and Xuhui (2010) methods concerned to different Checkland’s class can be interchangeable. So, techniques used in engineered system could be useful in natural systems. Mathematical models are formal statements or equations to express the relationship between system inputs, outputs and operations, to accomplish hydrological simulations. Then, mathematical models, as reasonable representations of the system, provide solution approaches to a wide spectrum of cases. For example, in streamflow or spring forecasting, recovering values in missing data, quantifying the impact over land changes, planning, designing and managing water resources, better understanding the hydrological processes, among others. The term Systems Engineering is a generic term that describes the application of structured engineering methodologies to the design complex systems which require high dynamic performance. Methods of system engineering are interdisciplinary tools to emulate complex systems; therefore multidisciplinary applications are required. These methods involve elements of mathematics, physics and computing as well as techniques of analysis and control, to provide solutions in many different fields; e.g. electronics, industrial mechanics, distributed computing, power, communication networks, manufacturing, logistics, artificial vision, robotics, transportation, chemical 27 1.3 Basic concepts INTRODUCTION processes, medical and biological systems, environmental systems, and bioprocesses. Fig. 1.2: Checkland’s taxonomy (Zexian and Xuhui, 2010): natural system, designed physical system (engineered system), designed abstract system and human activity system. Climate and karst (or glacier) are coupled systems in the context of system theory. Natural systems can leverage the tools from systems engineering Since general theory of systems attempts an integrative methodology for the treatment of scientific problems (von Bertalanffy, 1950), natural systems can gain the transferability of models from different scientific continents. The study of natural systems in this thesis leverages the systems engineering tools following the system identification techniques (Fig. 1.2). So, karst and glacier systems, from a data-driven analysis, are going to be characterized and identified by engineering methods essentially from a nonlinearity approach. From the perspective of systems engineering, the karst or glacier system can be assumed as a black box device with behavior ruled under physical laws, which are not fully known, but observations of the inputs and outputs can be useful in formulating a specific idea of the system dynamics. 28 INTRODUCTION 1.4 Hypothesis 1.4 Hypothesis System identification researchers should make decisions based on the following perspectives to cope with characterizing or predicting the behavior of a system based on recorded data: •How can they analyze the signals to obtain as much information as possible about the system? •How can they best use the information in the observed data to calculate a model with the same properties and behavior of the real system? •How can they know if the model is any good and how can I rely on it for simulation, design or prediction purposes? In the case of karst systems and glacier systems, these three questions are kept open to compelling answers, though already there are many and excellent approaches. As a starting point of this thesis to formulate the targets, and carry out a formal study applying the scientific method at each stage, I had to ask myself some questions on observed facts to pose well the problem, such as: •What is the parametric interdependence between input and output time series? •How can we unveil features in the signals of an aquifer that otherwise would remain hidden by methods heretofore known? •Why, in the same recharge area, with the same precipitation, does an aquifer present different flow regime than another? •Can I predict the discharge effects next time it rains in the recharge area of Ciano Polje (San Julian de Banzo)? •Does the complexity of nonlinear models in its theoretical structure and its computational implementation compensate improvements over linear models? •How can be quantified the interaction between climate and hydrological cycle by correlation analysis? •To what extent can a linear model emulate the dynamic behavior of a glacier? •Can be calculated the bounds of the active period in a glacier? •Which could be the suitable nonlinear structure for the annual cycle evolution of a glacier? •How can the glacier dynamics be classified by the features of glacier in the active periods? •Would be possible to achieve a generalization of the performance of a particular glacier in term of global model? 29 1.6 Precedent context INTRODUCTION •To what degree a glacier can play the rol of change climate sensor? 1.5 Aims of the thesis The primary objective of this thesis can be stated as: Applications of linear and nonlinear methods of system engineering to natural systems, namely karst aquifer and antarctic glacier, in order to reveal the relationship between variables through the spectral analysis, searching the best model for simulation and prediction, and contributing with advanced solutions that help optimize the management of a karst aquifer and characterize the hydrological dynamics of a glacier. The scope of this study has been the design and verification of models based on some linear techniques, with the support of spectral analysis based on wavelets techniques. Likewise, a set of selected nonlinear model structures has been used to the identification processes, focused on two natural scenarios previously chosen. One of them is the karst aquifer of Fuenmayor (San Julian de Banzo, Huesca, Spain), for the treatment of data collected from the years 2002 to 2005, and another is on Collins Glacier (King George Island, Antarctica) with data in the period 2001-2011. 1.6 Precedent context The Group of Technologies in Hostile Environments (GTE) from the University of Zaragoza (Spain) was founded in 1997, and its work is mainly oriented towards natural and/or hostile environments, such as high mountains, polar zones, canyons and underground environments (caves, mines, tunnels). The GTE, a multidisciplinary group integrated by researchers and lecturers, provides the proper background to develop a thesis like this. One of its research lines deals with the study of karst springs, mainly through a new system identification techniques, with results gathered in papers that come from over 10 years ago. Glackma Foundation promotes scientific research in the polar regions. Its researchers have been visiting both poles since 1985, almost each year. They have registered long time series related to glaciers discharge in Arctic and Antarctic glaciers, to study of the evolution of global warming, using glaciers as natural sensors. Glackma Foundation has provided the data and the fundamentals to carry out the part of this work about glaciers. 30 INTRODUCTION 1.7 Outline of the thesis 1.7 Outline of the thesis After exposing objectives, justification and fundamental assumptions of this work in the current chapter, this thesis is organized according to this scheme: Chapter 2 is a description of the main tools from system engineering that will be applied in the cases under study. It presents the theoretical basis to understand the fundamental of system identification techniques. It starts with classical analysis techniques and an overview of some existing methods for linear and nonlinear systems, with an introduction to an important block structured nonlinear system. The chapter also deals with specifically wavelets techniques, showing the powerful advantages of these mathematical tools that are going to be used extensively in the experimental cases considered. The last sections explain parametric and non parametric identification and finally the nonlinear identification is presented. Chapter 3 is devoted to the presentation of two natural systems to understand details of the conceptual models, identify the characteristics of the environment that condition the model design, and determine the scope of the problem to be worked out. The first review is focused on the karst system scenario, with a description of the karst hydrology and the input/output variables of the problem. The second one describes the antarctic environment, especially the melting processes in the glaciers. In order to understand weather and discharge correlation models, some fundamentals should be set. The chapter highlights the natural system features in both systems, and reviews the precedent models in the literature, in order to address the analysis and identification in next chapters. Chapter 4 begins with the description of the karst aquifer of Fuenmayor site. It follows analysis and system identification results after applying Linear time-invariant system theory, then nonlinear methods based on wavelet transform are applied to rainfall and discharge of aquifer, following some chosen techniques from Chapter 2 to get a system nonlinear model. Chapter 5 is the study of glacier Collins in Antarctica, drawing on fundamental techniques that have already been depicted in Chapter 2. It details some levels of coherence to work out the suitable seasonal effective for several years. Procedures based on linear methods have been tested, along with nonlinear tools, specially block models and analysis based on wavelet transform. Chapter 6 summarizes the most important conclusions of the results that have been reached throughout this thesis. It highlights the mostly research contributions and possible extensions of the research. Also, it discusses suggestions for further developments and improvements. 31 1.7 Outline of the thesis INTRODUCTION Finally, complementary reports are contained in appendices. •Appendix A contains the main abstracts of papers published in international journals and congresses written and presented by the author of this thesis or in collaboration with other researchers. •Appendix B contains examples of analysis results in the form of illustrations, which because of their complementary character, or over sizing, are not included in the text of the corresponding chapters. •Appendix C is a glossary about main terms defined according to the meaning given in the thesis context. 32 Chapter 2 System Identification Techniques An intellect knowing at any given instant of time, all forces acting in nature, as well as the momentary positions of all things of which the universe consist, would be able to comprehend the motions of the largest bodies of the world and those of the smallest atoms in one single formula, provided it were sufficiently powerful to subject all the data to analysis. Pierre Simon Laplace (1749-1827) 2.1 Introduction to System Identification Modeling is the abstraction of a real process to characterize its behavior. The scientific modelling aims to enhance the investigation of the phenomena in order to reveal and better understand cause-effect relationships (Williams, 1988). The model definition given by Eykhoff (1974) introduces the concept of “essential aspects”: “... [model] is a simplified representation of the essential aspects of an existing system (or a system to be constructed), which presents the knowledge of the system in a usable form”. The set of processes in a system determines the behavior of the system. Every process is determined by its physical and chemical properties, which are not always easily known. A model tries to emulate the “essential aspects” of the system behavior, simplified by choosing the most significant properties. So, modeling techniques can be classified as: •a priori modeling, white-box or morphological modelling, by making simple experiments to inquire the physical or chemical laws involved. •a posteriori modeling or black-box modelling, by building a model based only on data (data-driven) without having previous knowledge of the system. The model 33 2.1 Introduction to System Identification TECHNIQUES describes how the outputs depend on the inputs, not how the system actually is, and characterize the system dynamics (delays, speed, oscillations, and others), though the physical interpretation of the results is not straightforward. •the grey box modeling is an intermediate technique when peculiarities of internal laws are not entirely known, so it is based on both insight into the system and on experimental data analysis. System identification tries to estimate a black or grey model of a dynamic system based on observing input-output from experimental data. Zadeh (1962) defined system identification as: “... the determination on the basis of input and output, of a system (model) within a specified class systems (models), to which the system under test is equivalent (in terms of a criterion)”. The availability and reliability of the design techniques of system identification did expand the application fields beyond the scope of industrial applications. So that, there are models from system identification applied in other diverse fields, for example, economy, environment, biology, psychology, biomedical research, hydrology, and glaciology. Since the identification problem requires the model structures, a validation criterion and an aim (Ljung and Glad, 1994). Criteria and models will be exposed along this chapter. Some examples about identification aims could be listed here: •To design control strategies for a particular system (e.g., in optimizing an electrical microgrid operation). •To analyze the properties of the system (e.g., quantity rates in a medication reaction). •To forecast the evolution of the system (e.g. the future climate prediction according a IPCC downscaling model) •To identify hidden factors influencing in a system (e.g., sun spot in the karst spring). •To improve the internal knowledge of the system (e.g., the delay in the aquifer discharge respect to precipitation events). •To identify the interaction between coupled system (e.g., climate and glaciers). The objective of this chapter This chapter surveys the main approaches for identification and analysis of systems and their theoretical basis, in order to set the methodological scene for the study of the experimental cases described in later chapters, specially highlighting the procedures that mostly have been implemented in natural systems. 34 TECHNIQUES 2.1 Introduction to System Identification Outline of this chapter. The first section (2.2) of this chapter deals with problems in the acquisition of data from sensors as the sampling chosen (2.2.1) and treatment of outliers (2.2.2) to reach a sufficient time series quality for the subsequent treatment of the information. Sec. 2.3 is a review of the classical time series as the first approach to understand the underlying mechanism in the system. Two scenery are presented: time domain analysis (2.3.1) with autocorrelation structures, and frequency domain analysis (2.3.2) with an explanation of concepts like Fourier transform and frequency spectrum. Although wavelets techniques are included in the frequency domain issues, a special section is dedicated to wavelet techniques (2.4), because of its prominently application in the practical cases. The section contains a description of mother wavelets, detailed clarification of wavelet transform (2.4.2), process of estimation of wavelet power spectrum (2.4.3), and some discussions in using these techniques (2.4.4). A review of model structures (2.5) is necessary in the identification process to select suitable, identifiable model structure. Some linear time-invariant models are described in Sec. 2.5.1. Nonlinear models (2.5.3) go through Volterra series (2.5.3.1) and Hammerstein-Wiener models (2.5.3.2). Sec. 2.6 deals with several identification techniques, starting with the elucidation of the posed problem (2.6.1) and a simple literature overview (2.6.2). Parametric identification (2.8) relies on a model previously defined by a set of parameters that must be calculated to accomplish a given quality criteria. Linear techniques are given in (2.8.1) with advises on selection and verification criteria of the models (2.8.2). Nonparametric identification methods are described in 2.7, starting with methods in the frequency domain (2.7.2), as classical spectral analysis (2.7.2.1), going on with wavelet cross spectrum (2.7.2.2), and finishing with wavelet coherence fundamentals (2.7.2.3). In the time domain section (2.7.1) there are methods for system identification, as cross correlation (2.7.1.1) and impulse response (2.7.1.2). A final section provides an overview of nonlinear identification techniques (2.9) that is a tour around nonlinear parametric identification (2.9.1) remarking main issues on Volterra identification (2.9.1.1), and Hammerstein-Wiener identification (2.9.1.2). Also, methods about nonlinear nonparametric identification are mentioned (2.9). Finally, in the conclusion section (2.10), there are important assertions related to this chapter. 35 2.2 Time series TECHNIQUES 2.2 Time series To capture critical information about the processes to be investigated, field data are achieved through the sensor network. Process variables should be sampled for a duration and sampling frequency enough to obtain those quality time series that the analysis requires. The sequence of observations on one variable y(t), t ∈T (T is the discrete times domain), is called time series. The observation are usually equally spaced and indexed by integers (t = 1,... ,n) where n indicates the number of observations. The main objective of time series analysis is to get mathematical inferences from the sample data obtained from field sensors. A lot of problems can be solved by time series as: prediction (e.g. future weather precipitation), identification or abnormal peaks (e.g. outbursts in a glacier discharge), trends (e.g. sea global warming in the last century), etc. Before performing the analysis of the time series, some issues should be faced about the form of achieving raw data, check the integrity and reliability of preparing procedures; i.e., the sampling method, the outliers detection, and the lost data recovery. 2.2.1 Sampling period In many monitoring applications, there are some problem related to sensor devices to sample environment variables, mainly concerning with computing data as the memory size, processing capability and power supply. Sampling is the process by which continuous time signals –such as air temperature or water levels– are turned into discrete time signals. About batteries, if the sampling frequency is set too high, then the energy consumption would be so high that the sensor battery would be depleted too soon. To avoid this problem, the sampling frequency can be reduced, but this is not always possible. The Nyquist-Shannon sampling theorem states that if a function y(t) contains no frequencies higher than ωN, it is completely determinable by a sampling process of frequency 2ω. So, the sampling frequency ωSshould be bounded according the Eq. 2.1: 2ω=ωN<ωS<ωC(2.1) where ωNis the Nyquist frequency and ωCis the critical frequency for the sensor battery duration. On this matter, Alippi et al. (2009) presented an adaptive sampling algorithm for effective energy management in wireless sensor networks. Some of these sensors require computing capability to integrate a distributed artificial intelligence which offers a wide range of possibilities for the operation, automation and control of different systems (Hernández et al., 2012). 36 TECHNIQUES 2.4 Wavelet transform techniques Table 2.2: Requirements of a wavelet function ψ∈L2(Rd)(ˆ ψis the Fourier transform of the wavelet function ψ) Description Condition a The average value of the wavelet in the time domain should be zero ∞ R −∞ ψ(t)dt =0 b The function must have finite energy ∞ R −∞|ψ(t)|2dt =1 c Admissibility . The inverse wavelet transform only exists for 0 < Cψ<∞. This means that the analyzed signal can be reconstructed without loss of information. The constant Cψis called the admissibility constant. Cψ=2π ∞ R −∞ |ˆ ψ(ω)|2 |ω|dω The function ψ(t)is called «mother wavelet» or «basic wavelet» while the dilated and translated functions derived from the «mother wavelet» are called «daughter wavelets» or simply «wavelets» (Fig. 2.1). These daughter wavelets have the same shape as their mother wavelet. Their amplitude should rapidly decay away from the center of the wave in both time and frequent domains. The functional relationship between daughter ψs,τ(t) and mother ψ(t), in the scale s and displacement τ, is expressed as (Eq. 2.8): ψs,τ(t) = 1 √sψ(t)t−τ s(2.8) where s, τare real and s>0. Wavelets expressed by (2.8) include an energy normalizatation s−1/2which keeps the energy of the daughter wavelets the same as the energy of their mother. Morlet wavelet is the most popular complex wavelet used in practice, which mother wavelet is defined as Eq. 2.9 and represented in Fig. 2.1 ψ(t) = 1 4 √πejωt−e−ω2 2e−t2 2(2.9) where ωis the central frequency of the mother wavelet. Note that the term e−ω2 2is used for correcting the non-zero mean of the complex sinusoid, and it can be negligible for 5<ω For practical purposes, for large ω, e.g. ω>5, Eq. 2.9 can be simplified (FoufoulaGergiou and Kumar, 1994) by taking a complex cos wave modulated by a Gaussian envelope (Eq. 2.10): ψ(t) = 1 4 √πejωte−t2 2(2.10) Morlet wavelet provides a better energy localizing and higher frequency resolution, although the frequency-coordinate window shifts along frequency axis as scaling. 43 2.4 Wavelet transform techniques TECHNIQUES Fig. 2.1: Morlet’s «wavelet daughters» ψs,τ(t), that are dilated and translated functions derived from the «mother wavelet» ψ(t), according with Eq.2.8 2.4.2 Wavelet transforms The basic aim of wavelet analysis is both to determine the frequency (or scale) content of a signal and to assess and determine the temporal variation of this frequency content (Heil, 1989). Although wavelet analysis covers a wide range of methods and applications, fundamental operations are wavelet transforms, which are appropriate for the many natural phenomena that have the property that high frequency events happen for short durations. 2.4.2.1 Continuous wavelet transform (CWT) A wavelet transform correlates the signal with a family of waveforms ψs,τor wavelets –it is also called time-frequency atoms by Mallat (1998) and kernel by other authors– that meet the conditions of Table 2.2. The corresponding continuous wavelet time-frequency 44 TECHNIQUES 2.4 Wavelet transform techniques transform of f∈L2(R)is expressed by (Eq. 2.11): Wψ,s,τ{f(t)}=Z∞ −∞ f(t)·φ∗ s,τ(t)dt =hf,φs,τi(2.11) where Wψ,s,τ{f(t)}is the CWT of f(t) with basis function family ψs,τ. The wavelet coefficients represent a measure of similarity in the frequency content between a signal and a chosen wavelet function. These coefficients are computed as a convolution of the signal and the scaled wavelet function, which can be interpreted as a dilated band-pass filter because of its band-pass like spectrum. 2.4.2.2 Continuous wavelet transform with discrete coefficients (CWTDC) Signals are usually band-limited, which is equivalent to having finite energy, and therefore just a constrained interval of scales is useful. However, the continuous wavelet transform produces redundant information when capturing all the characteristics of the signal. The Discrete Wavelet Transform (CWTDC) has been created to minimize the redundancies produced by CWT. It is possible to compute the wavelet transform for just a proper selection of values of the frequency and time parameters and still not loose any information as recovering the original time series from its transform. The continuous wavelet transform with discrete coefficients is very similar to the continuous wavelet transform, but the parameters s and τare fixed to the power of 2, following these expressions: τ=2−jn,s=2−j,for m ≥0 and n ∈(−∞,∞)(2.12) The functional relationship between daughter ψk,j(t)and mother ψ(t), in the scale k and displacement j, is expressed as : ψk,j(t) = 2j/2ψ(2j−k)(2.13) So, Eq. 2.13 is a special case of Eq. 2.8. From the continuous wavelet function (Eq. 2.11) and the new values of s and τfrom Eq. 2.12, CWTDC takes the form of Eq. 2.14. Wψ,k,j{f(t)}=√2j ∞ Z −∞ f(t)·ψ∗ k,j(2jt−k)dt (2.14) 45 2.4 Wavelet transform techniques TECHNIQUES 2.4.2.3 Discrete wavelet transform (DWT) Although the output of the continuous wavelet transform contains discrete coefficients, its implementation could be hard, since the input signal is continuous. Discrete wavelet transform is the alternative. If f∈L2(R)and h∈L1(R), the convolution between the two signals: g(x)=(f⊗h)(x) = R∞ −∞f(t)h(x−t)dt. The continuous wavelet transform of a f(t) signal, given at Eq.2.11, yields a infinite set of wavelet coefficients. In a discrete form, where a time series uis a discrete sequence values of (u1,u2,...,un)separated in time by a constant time interval δt, the expression for the wavelet coefficient is given at time index j and scale s in the expression 2.15. Wψ,s,j(un) = N−1 ∑ n=0 unψ∗(n−j)δt s(2.15) where unis the discrete sequence, N denotes the length of the studied time series, ψ∗is the wavelet complex conjugated, and δtdenotes the sampling period. The algorithm to calculate the wavelet transform from Eq. 2.15, is a loop of convolutions performed N times for each scale, where N is the number of elements in the time series. The convolution theorem states that the Fourier transform of the convolution of two functions is the product of Fourier transforms of each function F{u(t)⊗y(t)}=F{u(t)}⊗ F{y(t)}. Performing the Fourier transform in both sides of the Eq. 2.15, the inverse of wavelet transform can be expressed as Eq. 2.16. Applying the Fourier inverse transform to FWψ,s,j(un), and following the work of Grinsted et al. (2004) , an efficient algorithm has been created to compute the discrete wavelet transform for the natural time series of this dissertation. This DWT can be computed by the inverse Fourier transform, as the Eq. 2.16 states. FWψ,s,j(un)= N−1 ∑ k=0 ˆuk·ˆ ψ∗(sωk)eiωknδt(2.16) where k = 0 ... N -1 is the frequency index, the Fourier transform of the function ψ(s,t)is ˆ ψ(sωk). The angular frequency ωkadmits the following values: ωk=(2πk)/(Nδt)for k≤N/2 ωk=−(2πk)/(Nδt)for k>N/2 2.4.2.4 Ouliers detection in the frequency domain Since outlier is an observation with characteristics of high-frequency phenomena, then wavelet technique is an excellent tool as outlier detector because of good location of 46 TECHNIQUES 2.4 Wavelet transform techniques frequencies. Most of energy and information of the data are usually concentrated in the first few coefficients. The outliers, as the noise, are on high frequency bands and reside in high-order coefficients. Therefore, the true data and outliers can be separated in the wavelet space (frequency domain). In a discrete wavelet analysis, a signal u(t) can be represented by a decomposition of the signal into approximations Ajand detailed coefficients Dj. This is accomplished using shifted and scaled versions of the original (mother) wavelet as given in Eq. 2.15. In practice, the wavelet coefficients are computed efficiently using the pyramid algorithm, introduced in the context of multiresolution analysis by Mallat (1998), that is based on a pair of high and low pass filters. The DWT in Eq. 2.15 produces a matrix of coefficients [ci,k]. Then, the function can be represented by: u(t) = j ∑ i_−∞ ∑ k ci,kψi,k(t) = j−1 ∑ i_−∞ ∑ k Di,kψi,k(t)+∑ j Aj,kψj,k(t)(2.17) Aiis the approximation component set of the time series, and represents the lowfrequency content. This low pass filter is like to continuously calculate a moving average of weighted data. Diis the detailed component set of the time series, and represents the high-frequency content. This high pass filter consists on a moving difference of the data. These two sets of wavelet coefficients facilitate the recursive form of the pyramid algorithm (Mallat, 1998). Struzik and Siebes (2000) propose a methodology capable of determining the statistical nature of the non-stationary process. The method checks the internal consistency of the scaling behavior of the process within the paradigm of the multifractal spectrum. Deviation from the expected spectrum is interpreted as the potential presence of outliers. Chinarro et al. (2011) proposed a wavelet-rosner test and applied it to the time series of a karst aquifer. The method also stems from the wavelet based multiresolution analysis. After checking the Aiand Dicoefficients (Eq. 2.17) to detect outliers, Rosner’s test is applied to remove or replace the abnormal values. A brief description of method in steps is: 1. Firstly, to avoid null values, a temporary shrunk series has been created with all elements from the raw data except null elements. 2. Transform the time series to the wavelet domain. DWT decomposes the time series and computes the approximation coefficients vector Aiand detailed coefficients vector Di, at level 1, following (2.17) 3. Apply the Rosner test on the Aicoefficients to get outliers in the frequency domain. On this step, another method to remove outliers can be used, but Rosner’s test has been weel tested in hydrological series. 47 2.4 Wavelet transform techniques TECHNIQUES 4. Eliminates all outliers from the Aiand from analogous index in Di, then two shrunk vector Arand Drare created. 5. To restore the time series, use Arand Drto compute the inverse wavelet transform. 2.4.3 Wavelet power spectrum (WPS) The wavelet power spectrum helps to estimate the repartition of energy in the signal to determine the concentration of a signal in singular instants and frequencies. Temporal variation in the distribution of energy across scales is one of the most usual applications of the wavelet transform. The Wiener-Khinchin theorem states that energy spectral density of a function is the Fourier transform of the corresponding autocorrelation sequence (Ricker, 2003). Analogously, wavelet power spectrum Pf(s,τ)is defined as autocorrelation function of the wavelet transformation (Wψ,s,τ)of u(t) and describes the power of the signal u(t) at a certain time τon a scale s: Pψ,s,τ{u(t)}=Wψ,s,τ{u(t)}∗W∗ ψ,s,τ{u(t)}=Wψ,s,τ{u(t)}2(2.18) Because the wavelet function ψ(τ)is in general complex, the wavelet transform Wψ,s,τ is also complex, with a real part, ℜ(Wψ,s,τ), an imaginary part, ℑ(Wψ,s,τ), an amplitude Wψ,s,τ, and a phase, ℜ(Wψ,s,τ)/ℑ(Wψ,s,τ). Analogously, WPS can be expressed in the same components. Wavelets add a new dimension in the spectral analysis to work simultaneously with time and frequency. Wavelet power spectrum is, in fact, a three-dimensional depiction, with time on the x-axis, frequency or scale on the y-axis, and the z-axis is to render the power magnitude at a particular time and frequency. This is a suitable tool for the spectral analysis of a non-linear system. 2.4.4 Wavelet transformation caveats 2.4.4.1 Cone of influence (COI) When applying the CWT to a finite length time series, the scalogram inevitably suffers from border distortions. The cause is that the values of the transform at the ends of the series cannot be accurately calculated because the transform calculus takes values 48 TECHNIQUES 2.4 Wavelet transform techniques outside the series range. These hedge effects also increase with s in a rate that depends on the mother function. The region in which the transform suffers from these edge effects is called the cone of influence (COI) and should be marked to take care in interpreting the belonged values (Mallat and Zhong, 1992). One solution is to pad the end of the time series with zeroes before applying the wavelet transform and then remove them afterward. The padding is an extension of time series should be sufficient to spread out the time series to the next power of two. The zero padding reduces the variance, but introduces discontinuities at the endpoints and decreases the amplitude near the edges as going to larger scales. (Torrence and Compo, 1998). As wavelet coefficients at COI suffer the same input discontinuity, the solution may be to rescale the remaining wavelet with choosing the e-folding time, i.e. as the distance at which the wavelet power drops by a factor e−2. Larger e-folding time implies more expansion of the wavelet spectrum, The e-folding time is a measure of the wavelet width, relative to the wavelet scale s and ensures that the edge effects are negligible above a threshold for a given signal u(t) (Grinsted et al., 2004). 2.4.4.2 Choosing wavelet function One singular characteristic of wavelet analysis is the arbitrary choice of the wavelet function. The below list is based on factors given by Torrence and Compo (1998), in order to select a suitable wavelet to get best performance, and has been completed with other considerations by the author of the thesis: a). Discrete or continuous. DWT provides a more compact representation of data. So, DWT is rather suited for image processing, signal coding, noise reduction and computer vision. Nevertheless, CWT (also CWTDC) is a transformation that provides a high redundancy of data, and is suitable for time series analysis and feature extraction purposes. b). Orthogonal or nonorthogonal. The use of an orthogonal basis implies the use of DWT while a non-orthogonal wavelet function can be used with either the discrete or the continuous wavelet transform. Orthogonal wavelet functions have a zero correlation each other while non-orthogonal wavelets have a nonzero correlation. Using an orthogonal wavelet, the signal can be transformed to the frequency domain and then return to the time domain with a negligible loss of information. Orthogonal wavelet analysis is useful for signal processing because it gives the most compact representation of the signal. Non-orthogonal wavelets tend to surplus energy, because of overlapping, and require a normalization to optimize the acquisition of information. They are useful for time series analysis. 49 2.4 Wavelet transform techniques TECHNIQUES c). Complex or real. A complex wavelet function will return information about both amplitude and phase and is better adapted for capturing oscillatory behavior. A real wavelet function returns only power, but is useful in location of peak frequency. d). Width. A wide wavelet function will give good frequency resolution and a loss of time resolution while a narrow wavelet function will yield good time resolution but a poor frequency resolutions. e). Shape. The wavelet function should reflect the type of features to be presented. For time series with sharp jumps or steps, a box-like function would be better such as Harr’s wavelet proposed by Haar (1910). Nevertheless, for smoothly varying time series, the recommendable wavelet is a smooth function such as a damped cosine. If a wavelet power spectra has to be performed, then ”the choice of wavelet function is not critical, and any function will give the same qualitative results as another” (Torrence and Compo, 1998). f). Choice of scales. In orthogonal wavelet analysis, the set of scales sis limited (Farge, 1992). In non-orthogonal wavelet analysis, an arbitrary set of scales can be built up for a more complete plotting. To analyze natural systems, this dissertation has chosen the Morlet wavelet in most cases, because of five interesting properties: •The peak frequency, the energy frequency and the central instantaneous frequency of the Morlet wavelet are all equal facilitating the conversion from scales to frequencies. •Heisenberg Box area has a reduced size with this wavelet, i.e. the uncertainty reaches a minimum value. Then Morlet wavelet has an optimal joint time-frequency concentration. •The time radius and the frequency radius are equal; therefore, this wavelet represents the best compromise between time and frequency concentration. •Finally, Morlet is a wavelet transformation that yields complex coefficients, with information on both the amplitude and phase. This facilitates the study of gaps and delays between two time series. 2.4.4.3 Advantages of wavelet transform Some features can be observed in the application of wavelet transform and according to Strang (1993) and Perrier et al. (1995). Nonlinearity. Analysis with Fourier transform is not completely successful in all types of problems. Exceptions are nonlinear systems, with very brief signals or sudden changes, 50 TECHNIQUES 2.4 Wavelet transform techniques as the typical time series given in a karst system and glacier discharge; hence, the study of their behavior should be carried out with different tools. Stationarity. Most traditional mathematical methods that examine periodicities in the frequency domain, such as Fourier analysis, have implicitly assumed that the underlying processes are stationary. The wavelet transform is suitable for the analysis of nonstationary signals because it provides a better time and frequency localization properties, expanding time series into time frequency space, such that the intermittent periodicities can better be localized (Perrier et al., 1995). Global properties. A Fourier transform hides information about time. It proclaims unequivocally how much of each frequency a signal contains, but is unknowable about when these frequencies were emitted. Therefore, in Fourier transform, any instant of a signal is similar to any other, even if the signal is as complex as public clap echoes in a theater, or changes as radically as the runoff through a river after a severe rainstorm. For an application, wavelets only capture the local time-dependent properties of data; whereas Fourier transforms, due to space-filling nature of the trigonometric functions, can only capture global properties (Popivanov and Miller., 2001). Computational efficiency Using the big Onotation (Becher et al., 2012), the computational complexity of the discrete Fourier transform is O(n2), where n is a number of time samples. Fourier transform is overtaken by the Fast Fourier transform (FFT) with a complexity O(nlogn)which takes less steps to solve an instance of the same problem. However, this is still under the complexity function O(kn)for discrete wavelet transform in decomposition and reconstruction processes (Fig. 2.2). The transform upshot with wavelets can be implemented in a computer by a quicker and more efficient algorithm. Fig. 2.2: Big-O complexity for Fourier Transform, FFT and DWT 51 2.5 Models TECHNIQUES 2.5 Models 2.5.1 Linear time-invariant (LTI) models Most of the properties of linear time-invariant (LTI) systems are due to the fact that the system can be represented by linear differential (or difference) equations. Such properties include impulse response, convolution, duality, stability, scaling, etc. The properties of linear, time-invariant system should not in general apply to nonlinear systems. Nevertheless, LTI could be a first approach to identify the non-linear system. The effect of any invariant linear system (LTI) on an arbitrary input signal is obtained by convolution of the input signal with the system’s impulse response function. In a LTI system, the output of the system y(t) for an input x(t) can be obtained by the convolution integral: y(t) = g(t)?x(t)Zt 0 g(t−τ)x(τ)dτ(2.19) where g(t) is the impulse response of the system. That is, g(t) is the output of the system with an input x(t) = δ(t), where δ(t)is the Dirac delta. The impulse response completely characterizes the dynamic behavior of the system. Applying the Laplace transform to the convolution integral (Eq. 2.19) we obtain Eq. 2.21 : L[y(t)] = L[g(t)∗x(t)] = L[g(t)]L[x(t)] (2.20) or in simple expression: Y(s) = G(s)X(s)(2.21) where Y (s), G(s) and X(s) are the Laplace transforms of y(t), g(t) and x(t) respectively. ATransfer Function (TF) is the mathematical representation of the relation between the input and output of a system. In a LTI system, TF can be expressed as the ratio of the Laplace transform of the output and the input, and corresponds to the Laplace transform of the impulse response G(s). G(s) = Y(s) X(s)(2.22) The transfer function of a system is rational fraction with numerator and denominator polynomials of the complex variable s: G(s) = bmsm+...+b0 ansn+...+a0 =N(s) D(s)(2.23) The roots of N(s) are called zeros of the system and roots of D(s) are called poles of the system. Poles and zeros are complex numbers that determine the dynamic behavior 52 TECHNIQUES 2.7 Non-parametric identification 2.6.2 The literature highlights. This section can not describe all historical processes on the system identification, only some highlights are mentioned from the engineering views A further purpose can be found in the work of Deistler (2002), with an excellent review of the history of system identification and time series analysis. Spectrum analysis of time series might have commenced in 1664, when Isaac Newton decomposed a light signal into frequency components as passing the light through a glass prism. In 1800, Herschel measured the average energy in various frequency bands of the spectrum by placing thermometers along each band. The first foundations on identification processes were set by mathematicians – Gauss (1809) and Fisher (1912)–, although with subsequent important contributions from engineering and statistics. Aström and Bohlin (1965) introduced the Maximum Likelihood framework, based on projection techniques in Euclidean space, which have been used extensively on the estimation of the parameters of difference –also known ARMA (AutoRegressive Moving Average) or ARMAX model (AutoRegressive Moving Average with eXogeneous inputs). Among contributions of Ljung and Glover (1981), there is one to clearly separate two independent concepts: the choice of a parametric model structure and the choice of an identification criterion. In system identification, there are two approaches: parametric and nonparametric identification. In the parametric identification problem, a mathematical structure is assumed to govern the system, and the identification processes is focused only on the determination of unknown parameters for this structure that optimize the representation of the system. Nonetheless, in a non-parametric identification the structure of these equations is also unknown. Nonparametric regression and spectral techniques correspond to this kind of techniques (Box and Jenkins, 1976; Ljung, 1999). 2.7 Non-parametric identification Nonparametric identification techniques provide a very effective and simple way of finding model structure in data sets without the imposition of a parametric one. Commonly, the initial process to carry out is the nonparametric identification, and then, if it were suitable, the parametric identification should be performed. The next sections review the nonparametric identification methods from time domain and frequency domain perspectives. 59 2.7 Non-parametric identification TECHNIQUES 2.7.1 Non-parametric identification in the time domain 2.7.1.1 Cross-correlation Cross-covariance is a non-parametric identification technique and is related with the impulse response g(t) (Eq. 2.19) of a system and thus with its behavior (Box et al., 1994; Ljung, 1999). Assuming that the signals have zero mean, if y?(t)and x?(t)are uncorrelated, the correlation between the input and the output is: y?(t) = ∞ ∑ k=0 g(kT)x?(t−kT)+v?(t)(2.32) where v?(t)is the noise in the system. The signals involved can be regarded as the realization of stochastic processes. We can define the following coefficients and functions: If v?(t)and x?(t)are uncorrelated, the (cross) covariance function between the input and the output is: Rxy(τ) = ∞ ∑ k=0 g(kT)Rx(τ−kT ) = g?(τ)?R? x(τ) That is the cross correlation is the convolution between the impulse response and the autocorrelation of the input. Thus, the impulse response can be estimated from the covariance (correlation if both signals have zero mean) if the input is a white noise. However, this is not a common case. For example, in hydrological systems, we have no control over the time series that always differs from the white noise. This problem is solved by the use of a whitening filter over the input and the output (Box et al., 1994; Ljung, 1999). 2.7.1.2 The impulse response When the system has a finite impulse response, the non-parametric identification is performed by an intermediate approach between non-parametric and parametric identifications, and corresponds with a representation of the system by a FIR structure. Two methods have been proposed to calculate the coefficients of the impulse response: 60 TECHNIQUES 2.7 Non-parametric identification Method 1 is based on the Wiener-Hopf summation equation (Labat et al., 1999c, 2000a). The estimator of cross-correlation between the input x(t) and output y(t) of the system for an infinite time series, is given by: RN xy(i) = M−1 ∑ k=0 g(k)Rx(i−k)withi=0..M−1 (2.33) Deploying the equation 2.33 in a matrix form (Labat et al., 2000a):    Rx(0)Rx(−1)... Rx(1−M) Rx(1)Rx(0)Rx(−1)... ... Rx(1)Rx(0)Rx(−1) Rx(M−1)... Rx(1)Rx(0)   ·   g(0) g(1) ... g(M−1)   =   Rxy(0) Rxy(1) ... Rxy(M−1)   (2.34) Thus, the impulse response can be calculated as: Rx·g=Rxy ⇔g=R−1 x·Rxy (2.35) Method 2 is based on a minimization process. The most commonly used are linear programming (Dreiss, 1982,9) and least squares (Dreiss, 1989; Padilla and Pulido-Bosch, 1990; Deni´ c-Juki´ c and Juki´ c, 2003). To identify the impulse response is to minimize the prediction error. Considering a FIR model of length M: Q?(t) = M−1 ∑ k=0 g(kT)x?(t−kT)+ε(t) = ˆ Q?(t)+ε(t)(2.36) The minimization of ε(t)can be done by several methods. In karst hydrology, least squares and linear programming over the slack variables has been used. The minimization of prediction error of a FIR model is a case included in the so called parametric identification methods (Ljung, 1999). However, the parametric identification considers also several model structures with infinite impulse response and several models for error behavior. 2.7.2 Non-parametric identification in frecuency domain 2.7.2.1 Spectral analysis Frequency response can be derived from Fourier transform of the impulse response signal. Provides information about the gain and phase of the system for different input frequencies. The cross spectrum can be calculated in the form: 61 2.7 Non-parametric identification TECHNIQUES φN xy(n) = M−1 ∑ k=−M+1 w(k)RN xy(k)e−jn 2π Nk(2.37) Let g? T(t)the impulse response of a sampled system. Then, as it has been stated, the output y?(t)from an input x?(t), follows: y?(t) = ∞ ∑ k=0 gT(kT)x?(t−kT)+v?(t)(2.38) where v?(t)is the noise in the system. If v?(t)and x?(t)are uncorrelated, the (cross) covariance function between the input and the output is: Rxy(τ) = ∞ ∑ k=0 g(kT)Rx(τ−kT ) = g?(τ)?R? x(τ) So, the cross correlation is the convolution between the impulse response and the autocorrelation of the input. Supposing the sequences of finite length N, applying the Fourier transform and the equations 2.7 and 2.37 we obtain: φN xy(ω) = GN(jω)φN x(ω)(2.39) thus, an estimate of the frequential transfer function can be obtained as GN(jω) = φN xy(ω) φN x(ω)(2.40) However, the result is a plot which is not recommendable directly for simulation purposes (Ljung and Glad, 1994). Another useful function derived from the cross spectrum is the Coherence Function, which is calculated by the following expression: γ2 xy(ω) = |φN xy(ω)|2 φN x(ω)φN y(ω)(2.41) It measures the linear correlation between the input and the output of the system at each frequency ω. Notice also that the Coherence Function is dimensionless and can be shown that 0≤γ2 xy(ω)≤1. 62 TECHNIQUES 2.7 Non-parametric identification 2.7.2.2 The cross wavelet spectrum (XWS) The bivariate extension of wavelet analysis is recommended when the system involves two time series, instead of only one, to assess time-varying spectral relations between two signals which are often non-stationary. Cross-wavelet transform (XWS), – also called wavelet cross-scalogram or coscalogram–, gives information on the dependence between two signals, u(ti) and y(ti), as a function of time, similar to cross-correlogram (or the cross-spectrogram). XWS, a bivariate extension of WPS, is commonly applied in earth sciences, e.g., in the experimental cases of this document, the air temperature and glacier discharge or rainfall and conductivity in an aquifer system. The cross wavelet spectrum, introduced by Hudgins et al. (1993) to study the atmosphere turbulence, reveals how regions in the time frequency space with large common power have a consistent relationship. This fact suggests a causality between both time series, and the expectation value of the correlation of two signals u(t) and y(t) is: Xψ,s,τ{u(t),y(t)}=W∗ ψ,s,τ{u(t)}∗Wψ,s,τ{y(t)}(2.42) Since the wavelet function ψ(τ)is in general complex, Xψ,s,τ{u(t),y(t)}is also complex. The cross wavelet power of two signals describes the covariance between these times series at each scale of frequency. Cross wavelet spectrum illustrates quantitatively the similarity of power between two times. It has already been applied to rainfall-runoff cross analysis by Labat et al. (2001), and was briefly discussed by Maraun and Kurths (2004). If complex wavelets are used, such as the Morlet wavelet, its squared absolute value |Xψ,s,τ{u(t),y(t)}|2or simply its absolute value |Xψ,s,τ{u(t),y(t)}|is often gotten for better plotting. The value of |Xψ,s,τ{u(t),y(t)}|2is large when |Pψ,s,τ{u(t)}|and |Pψ,s,τ{y(t)}| are big at the closeness of scales (frequencies) and around the same time, regardless of the local phase difference. When the phase information is required, the Eq. (2.42) should be expressed in terms of its module and phase angle: Xψ,s,τ{u(t),y(t)}=|u(s,τ)|e−iθu(s,τ)|y(s,τ)|eiθy(s,τ)=Xψ,s,τ{u(t),y(t)}eiθy(s,τ)−iθu(s,τ) (2.43) This means that the phase angle iθy(s,τ)−iθu(s,τ)reflects the phase difference by which y(t) leads u(t) at the given scale and time. Van Milligen et al. (1995) introduced delayed wavelet cross spectrum, a useful quantity to detect structures from two separated observation points. 63 2.7 Non-parametric identification TECHNIQUES In hydrology, wavelet cross-correlation should sometimes be preferred to classical cross-correlation, because this method provides new insights into the scale depending on the degree of correlation between two given signals (Onorato et al., 1997; Labat et al., 2000b). 2.7.2.3 Wavelet coherence spectrum (CWS) An extension to Fourier analysis, to allow for non-stationarity, is windowed Fourier analysis, but this overcomes the assumption of global stationarity within each segment. In the time-scale domain, cross-spectrum cannot be normalized locally assuming stationarity to have a value bounded by, for example, zero and one, because multiple points should be involved for some degree of smoothing. Coherence in signal processing consists of a measure of the correlation between two signals. Power Spectrum represents the power carried by each frequency in a signal. The similarity of two signals can be checked by estimating CWS. Two examples coherence application are: the wavelet local correlation coefficients introduced by Buresti and Lombardi (1999) to measure the phase coherence of the signals; and cross wavelet coherence function introduced by Sello and Bellazzini (2000) to assess the intensity coherence of turbulent signals. In practice, to calculate the coherence of two signals, we calculate the cross-power spectrum, because the coherence is the normalized measurement of the cross-power spectrum and can be calculated dividing cross-power spectrum by the squared root of the product of the spectra of the signals (Eq. 2.44). The expression of wavelet coherence is (Liu and Todini, 2002): Cψ,s,τ{u(t),y(t)}=Xψ,s,τ{u(t),y(t)} Pψ,s,τ{u(t)}·Pψ,s,τ{y(t)}1/2(2.44) The polar expression is a useful form of wavelet coherence with modulus ρψ,s,τand phase φψ,s,τ: Cψ,s,τ{u(t),y(t)}=ρψ,s,τ{u(t),y(t)}·eiφψ,s,τ(u(t),y(t)) (2.45) A value of 1 means a perfect coupling between u(t) and y(t) around time τon a scale s for the wavelet ψ. For a zero or negative value the variations of two signals are not correlated, and for positive value between 0 and 1 the variations are correlated in a certain degree. This value meaning is similar to a traditional correlation coefficient as defined by (Barrett (1964)), and it is useful to think of the wavelet coherence as a localized correlation 64 TECHNIQUES 2.8 Parametric identification coefficient in time frequency space. Due to these properties, wavelet coherence is an increasingly popular method in analyzing hydrological correlated events. The fact a signal is correlated with other, not only means that some energy in a frequency is present in both signals, but that plus the frequency which is present in both signals is also related by phase. The magnitude squared coherence is also used, which is the squared value of the cross-power spectrum divided by the product of the power of the spectra of both signals. The measure of coupling average in the scale (frequency) domain between input and output signals u(t) and y(t) would provide information about the relationship of the signals. This can be got by the wavelet cross spectrum coefficient or simply wavelet coherence spectrum (CWS), obtained from the normalized wavelet cross-spectrum (to have values between 0 and 1). 2.8 Parametric identification Parametric identification relies on a model previously defined by a set of parameters that must be calculated to accomplish a given quality criteria. E.g., the system characteristics can have a parametric representation through a polynomial of a finite and known degree. The model structure can be obtained by physical modeling (grey box) or it may be a standard one (black box). In the latter case, a set of generic standard structures must be taken into account (OE, FIR, ARX, ARMAX and BJ). 2.8.1 Linear parametric identification Parametric identification techniques depend mostly on Prediction-Error Methods (PEM) (Ljung, 1999). The output of system y?(t)can be expressed as Eq. 2.32. A more useful expression is based on the Ztransform: Y(z) = G(z)X(z)+W(z)(2.46) This expression can be rewritten as: Y(z) = G(z)X(z)+H(z)E(z) = N(z) D(z)X(z)+ A(z) B(z)E(z)(2.47) where E(z)is the transform of a white noise, ε(t).G(z), the transfer function of the system, and H(z), the stochastic behavior of noise. G(z)and H(z)are rational functions 65 2.8 Parametric identification TECHNIQUES whose numerator and denominator are polynomials of the zvariable. The relationship between both functions defines several model structures. Figure 2.6 shows the most common ones: AutoRegressive eXogeneous (ARX) model, AutoRegressive Moving Average eXogeneous (ARMAX) model, Box-Jenkins (BJ) model and Output Error (OE) models. The features and advantages of each structure have been studied in several works, see for example Ljung (1999). The ARX model D(z)Y(z) = N(z)X(z)+E(z) is the easiest to estimate since the corresponding estimation problem is of a linear regression type. The foremost disadvantage is that the disturbance model 1/N(z) comes along with the system’s poles. It is, therefore, easy to get an incorrect estimate of the system dynamics because the A (Eq. 2.47) polynomial can also include the disturbance properties. So, higher orders in A and B coefficients in the Eq. 2.47 may be required. If the signal-to-noise ratio is good, this disadvantage is less important. The ARMAX model D(z)Y(z) = N(z)X(z) + A(z)E(z) gives extra flexibility to handle disturbance modeling because of the C polynomial (Eq. 2.47)). For this reason, ARMAX is a widespread used model and performs well in many engineering applications. The OE model Y(z) = [N(z)/D(z)]u(z) + E(z) has the advantage that the system dynamics can be described separately and that no parameters are wasted on a disturbance model. If the system operates without feedback during the data collecting, a correct description of the transfer function G(z) = N(z)/D(z) can be obtained regardless of the nature of the disturbance. In the BJ model Y(z) = [N(z)/D(z)]u(z) + [A(z)/B(z)]E(z) the disturbance properties are modeled separately from the system dynamics. 2.8.2 Selection and verification criteria To get a model reliable the results and predictions inferred from model should be verified and validated. Model validation is carried out by comparing the model behavior with the system’s one and evaluating the difference. All models have a certain domain of validity. This may determine how exactly they are able to describe the system behavior. Efficiency criteria can be defined as mathematical measures of how well a model simulation fits the available observations (Beven, 2001). A number of different methods to set a criterion have been suggested in the literature, e.g. least squares, generalized least squares, maximum likelihood or instrumental variables. Krause et al. (2005) have studied the utility of several efficiency criteria in three examples using a simple observed streamflow hydrograph, declaring that: “The selection and use of a specific efficiency 66 TECHNIQUES 2.8 Parametric identification OE :Y(z) = N(z) D(z)X(z)+E(z) FIR : when D(z) = 1 ARX :Y(z) = N(z) D(z)X(z)+ 1 D(z)E(z) ARMAX :Y(z) = N(z) D(z)X(z)+ A(z) D(z)E(z)BJ :Y(z) = N(z) D(z)X(z)+ A(z) B(z)E(z) Fig. 2.6: System model structures. criteria and the interpretation of the results can be a challenge for even the most experienced hydrologist, since each criterion may place different emphasis on different types of simulated and observed behaviors”. Nash-Sutcliffe efficiency, coefficient of determination, and index of agreement, are frequently applied to verify hydrologic models. The efficiency value E proposed by Nash and Sutcliffe (1970) is defined as one minus the sum of the absolute squared differences between the predicted and observed values normalized by the variance of the observed values during the period (Eq. 2.48). E=1−∑N k=1[y(kT )−ˆy(kt)]2 ∑N k=1[y(kT )−y]2(2.48) where y(kT)is observed data in the sampling interval T (with k=0,1,2,...,N), ˆy(kt)is modeled output, and yis mean of observed data. Nash-Sutcliffe efficiencies can range from −∞to 1. An efficiency value of 1 (E = 1) corresponds to a perfect match of model output to the measured data. An efficiency value of 0 (E = 0) indicates that the model is as accurate as the mean of the observed data, whereas an efficiency less than zero (E < 0) occurs when the observed mean is a better predictor than the model or, in other words, when the residual variance (described by the numerator in the expression above), 67 2.9 Nonlinear identification TECHNIQUES is larger than the data variance (described by the denominator). Essentially, the closer the model efficiency value is to 1, the more accurate the model is. According to Legates and McCabe (1999), the disadvantage of the Nash-Sutcliffe efficiency is the fact that the differences between the observed and predicted values are calculated as squared values. As a result, larger values in a time series are strongly overestimated, whereas lower values are neglected. E.g., in the quantification of an aquifer discharge, this criterion could lead to an overestimation of the model performance during peak flows and an underestimation during low flow conditions. Nevertheless, Nash-Sutcliffe efficiency, expressed frequently as 0-1 coefficient (E), is the unique criterion used in the experimental cases of this thesis, in order to standardize the comparison between models. 2.9 Nonlinear identification It is difficult to establish a clear identification methodology of nonlinear systems, since analysis is usually more intricate than in the identification of linear models, because of the variety of nonlinear model structures and nonlinear behaviors. Donoho and Johnstone (1994), and Donoho (1995) introduced nonlinear wavelet estimators in nonparametric regression through thresholding, i.e., the term-by-term assessment of coefficients in the wavelet expansion. Only coefficients that exceed a predetermined threshold are taken in account. This produces the wavelet shrinkage. Bendat (1990) describes procedures to identify and analyze the properties of many types of nonlinear systems as Zero-Memory Nonlinear Systems and Parallel Nonlinear System, with analysis of Nonlinear System Input/Output Relationships. Zhang (1997) applied wavelet theory for nonlinear system identification, with a wavelet basis as a universal function approximator, with a neural network used to determine the resolution, and the translation coefficients of the wavelet. This nonparametric estimator named wavelet network has a neural network like structure that makes use of techniques of regressor selection completed with backpropagation procedure. This section is going to focus only on nonlinear parametric identification, and, inside this type, Volterra series and HammersteinWiener methods. 68 KARST AND GLACIAL HYDROLOGY 3.2 Karst background 3.2 Karst background Karst word is referred to the German expression that is stemmed from the Slovene word kras and the Italian carso, and means rocky terrain. Karst is a landform originated primarily by dissolution on limestone, with diverse shapes, blind valleys, caves and singular springs. The term was originally used to describe a landscape in the western part of Croatia and Slovenia, where this kind of terrain was firstly studied and formally defined. Currently, the name is internationally accepted to designate similar processes anywhere in the world. Over decades, geomorphology and hydrology of karst was separately studied as different issues. However, from 1980 both perspectives are regarded together as this karst definition explains: “Karst is a terrain with distinctive hydrology and landforms arising from the combination of high rock solubility and well developed solutional channel (secondary) porosity” (Ford and Williams, 1989b). An emerging approach to karst hydrogeology, gives a perspective more integrative and universal by encompassing the whole range of karst processes as this declaration: “An integrated mass-transfer system in soluble rocks with a permeability structure dominated by conduits dissolved from the rock and organized to facilitate the circulation of fluid” (Klimchouk, 2004). In order to understand the groundwater aspects, the next section starts with a brief introduction to the karst processes. 3.2.1 Karst processes. Cvijic (1893) is considered the pioneer in understanding the genesis of the karst cavities. His work provides a relation between the existence of karst landscapes and the groundwater flow, describing typical vertical development zones (dry, transition and saturated). He stated that the dissolution of the rocks was the key process for the formation of caves and sinkholes. Later, Ford and Ewers (1978) elucidates the karst formation process, though it is more detailed in the publication of Ford and Williams (2007). Karstification mainly consists of a weathering due to dissolution in water, whose aggressiveness is characterized by the action of other chemical processes (hydration, ionic substitution and oxidation-reduction), and physical phenomena (mass transfer and diffusion). Weathering is defined by Summerfield (1991) as “the adjustment of the chemical, mineralogical and physical properties of rocks in response to environmental conditions prevailing at the Earth’s surface”. Parameters, that in general determine the 75 3.2 Karst background KARST AND GLACIAL HYDROLOGY weathering process, are: the type of geological material, the atmospheric pressure, temperature, and contact of water with air. The weathering process is part of the carbon cycle, in which carbon is exchanged between the atmosphere, surface water, groundwater and carbonate minerals. Limestone, the most common sedimentary rock in karst formations, is mainly composed of the mineral calcite, the most stable polymorph of calcium carbonate (CaCO3). A carbonate rock with at least by a 90% of calcite is considered a pure carbonate rock (Jennings, 1985). The limestone is produced by chemical precipitation and the accumulation of organisms such as corals and gastropods. The dissolution of calcium carbonate by carbon dioxide in aqueous solution is the dominant reaction in karst processes, including speleogenesis. The reaction can be represented by: CO3Ca →CO2− 3+Ca2+Γ(3.1) Along with the water, CO2plays an important role in the karstification process. It comes generally from the atmosphere, though also CO2can be produced biologically in the soil. CaCO3+CO2+H2O→Ca2++2HCO2− 3(3.2) The reaction is reversible. The solution, containing the dissolved calcium bicarbonate, can lose carbon dioxide to the atmosphere and precipitate calcium carbonate. This process is responsible for the development of speleothems, and tufa or travertine at the surface (EPA, 2002). Precipitation of dissolved carbonate minerals is accompanied by the release of the carbon as CO2. The dissolution progresses toward deep underground areas with aggressive water filtration through multiple channels such as potholes, sinks, joints and chasms. In the vertical penetration, the water is away from contact with atmospheric CO2, then the water aggressiveness lessens, losing the ability to dissolve. 3.2.2 Karst classification The classification of karst is very important from the standpoint of hydrogeology. Cvijic (1893) provided one of the first classifications of karst according to the degree of karstification: Holokarst (complete karst) and Merokarst that defines an imperfect karst topography Ford and Williams (1989b) classify karst process according to morphology, in endokarst and exokarst, although the author admits that there may be mixed formations. 76 KARST AND GLACIAL HYDROLOGY 3.2 Karst background •Exokarst. All features that may be found on a surface karst landscape, ranging in size from tiny karren to large poljes, and from protruding forms to deep depressions, belong to the exokarst (EPA, 2002). Exokarst can be prominent forms as Karren or depressions as Doline (a concave depression in limestone) and Polje that is a large depression in a karst which follows the main geological structure trends. •Endokarst. Endokarst is the layer of the karst system beneath the surface. It is created by surface waters that slowly percolates dissolving the limestone, creating deep crevices, galleries and large caves, where water flow arises to constitute the underground runoff. 3.2.3 Karst hydrology An aquifer (from aqua = water and fero = carry) is all geological formation saturated in water that is capable of acquiring, maintaining and passing water through geological formations. Accordingly, the aquifer comprises in a broad sense, a recharging surface, a vertical flow of water through the vadose zone, a storage and transmission, and a discharge area that can be very reduced. There are three types of aquifers: detrictic, fissured and karst. The last is the typical of the limestone terrains. From the notion of system, Yevjevich (1959) firstly used the concept of the hydrological system. Mangin (1975), leading the Laboratoire souterraine of Moulis at the French Pyrenees, translated it to the concept of «karst system». He defined karst system as the entity where the water flow constitutes a drainage unit and introduced an important paradigm shift on the study of karst hydrology, linking together karst genesis, morphology, hydrology and hydrochemistry. In practice, many authors use karst system as synonymous of karst aquifer, and also in this work. Bakalowicz (2005) described the various types of recharge and functions in the vadose and phreatic zones of a karst system. Karst aquifer, sometimes visitable, can exist even when there are no recognizable karst landforms on the surface, and even when there are no known and accessible caves. The knowledge of each aquifer requires a study of recharge, internal constitution, water flowing through the rocks, the chemical and physical processes between water and rocks, and finally the process by which the aquifer is discharged, either naturally through springs and rivers, either by artificial procedures as pumped wells. 77 3.2 Karst background KARST AND GLACIAL HYDROLOGY 3.2.3.1 Karst conceptual scheme Since the first karst studies, there was a controversy on the existence of a saturated zone. The French school, headed by E. Martel, denied the existence of a saturated zone, based on their caving experience in the Pyrenean underground rivers (Martel, 1910). The Austrian school, from the study of Slovenian karst, recognized the existence of underground water zones (Cvijic, 1893). Along the time, the ideas of both schools were going mingling till to converge, according to Ford and Williams (1989b). Underground water, as such, is hardly seen until it ceases to be underground. For this reason, it is difficult to study and measure (Bakalowicz, 2005). The phenomenon of the progression of the water flow in a karst system is complicated and apparently not completely known. In the literature, most of the authors (Dooge, 1973; Mangin, 1975) agreed, though differing in the terms used, to divide the karst system into several parts or subsystems to tackle the model study. Essential components of a karst aquifer and system structure were stated by White (2003) and sketched in Fig. 3.1. Fig. 3.1: The essential components of the karst aquifer (White, 2003) Currently, most conceptual models distinguish three main zones –besides superficial recharge area– in the vertical direction, taking into account the essential hydrogeological characteristics described by Mangin (1975): 78 KARST AND GLACIAL HYDROLOGY 3.2 Karst background •Epikarst zone, with a relatively thickness that may vary significantly (15 to 30 meters may be a good generalization), is the portion of bedrock that extends from the base of the soil zone and is characterized by extreme fracturing and intense solution. It is separated from the phreatic zone by an inactive, relatively waterless interval of bedrock that is locally breached by vadose percolation. The epikarst develops a higher activity, due to the proximity to the surface, where water has a higher carbon dioxide concentration. Significant water storage and transport are known to occur in this zone (Klimchouk, 2004). •Transit zone. It is below the epikarst, with which it forms the so-called unsaturated zone and is more defined by vertical lines that provide a rapid progress infiltration reaching the water table or piezometric. •Saturated zone or phreatic zone comprises those parts of the karst in which all voids are filled with water. The saturated zone is formed by conduits of large dimensions, narrow cracks and pores in the rock. In many karst conduits are organized in a hierarchical network similar to a river system. This network finishes at the surgence area before commented. In the conceptual model proposed by Mangin (1975), the main conduit system transports infiltration waters towards a karst spring, but is poorly connected to large voids in the adjacent rocks, referred to as the "’annex-to-drain system"’. The particularity of Mangin’s conceptual model is that it associates the storage function of karst aquifers to the "’annex-to-drain system"’. The conceptual model proposed by (Drogue, 1974) assumes that the geometric configuration of karst conduit networks follows original rock fracture patterns with a size in the order of several hundred meters, separated by high-permeability low storage conduits. The conceptual model proposed by Király (1998) combines both Mangin and Drogue models. Although Királys’ model employs a hierarchical conduit network similar to the model of Mangin and takes the effect of the epikarst into account, it also comprises the hydraulic effect of nested discontinuity groups similar to those involved in the model of Drogue. Király demonstrated the scale effect to be a consequence of coexisting discontinuity groups of different scale associate it to the low-permeability matrix. During last decades of twenty century, an important number of studies were accomplished on the karst hydrology. The result was the identification of several types of karst aquifers. White (2007) proposed a conceptual scheme for classifying carbonate aquifers in terms of groundwater flow system and hydrogeologic setting, recognizing three principal flow type: diffuse, with slow laminar flow; free, with fast and turbulent flow in conduits; and confined, with generally slow flow in thin layers intercalated between impervious rocks. 79 3.2 Karst background KARST AND GLACIAL HYDROLOGY 3.2.4 Modeling approaches The complexity in modeling the karst system, already mentioned by Lastennet and Mudry (1997), is determined by the high capillarity in the water flow scattered through ducts of an internal structure that is not well known (Bakalowicz, 2005). Attempts to model karst aquifers using spatially distributed numerical models proved acceptable results on a bounded local scale for management purposes (Angelini, 1997; Larocque et al., 1998). Nevertheless, this type of numerical model requires quite large databases (permeability, porosity, transmissivity, piezometry, etc.) and computing complexity, specifically when the system is heterogeneous, as the most karst aquifers are. There are two fundamentally different types of modeling approaches for quantitatively characterizing the hydraulic behavior of karst hydrogeological systems (Martin and White, 2008). Global models consist of the mathematical analysis of spring discharge time series (hydrographs). According to this approach, karst systems can be considered as transducers that transform input signals (recharge) into output signals (discharge). As the acquisition of spring discharge data is relatively simple, these models have been already used since the beginning of the last century. However, global models do not take into account the spatial variations within the aquifer. Consequently, they cannot provide direct information concerning aquifer hydraulic parameter fields. Distributive models are suitable for the quantitative characterization of the spatial variations of hydrogeological phenomena. Distributive methods consist of subdividing the model domain into homogeneous sub-units, and calculating groundwater flow by applying flow equations derived from basic physical laws. 3.2.4.1 Input: Recharge In general the aquifer recharge is the volume of water entering the aquifer, that mostly can be estimated by the relative percentage of precipitation in the basin, which constitutes the system input signal. Recharge can be estimated by different methods. E.g., Andreo et al. (1996) develop a simple method to estimate the rate of carbonate aquifer recharge (expressed as a percentage of precipitation) by combining different variables (geological, geographical, morphological and soil), and to set the zonal distribution of recharge aquifers. 80 KARST AND GLACIAL HYDROLOGY 3.2 Karst background There are two basic ground-water recharge types in karst terranes: autogenic (concentrated or diffuse) and allogenic (Shuster and White, 1971). •Autogenic recharge is derived directly from precipitation over a closed depressions (i.e. gorges, poljes or canyons), and enters the rock at discrete points (sinkholes or large fractures). The autogenic recharge is called diffuse when there is a slow percolation through numerous small orifices and in a large area of limestone outcrop. •Allogenic recharge to karst aquifers occurs when the surface runoff is on large areas of insoluble rock or low permeability soils, so the water flows directly from the adjacent soluble carbonate bedrock (Palmer, 2002). Effective rain. Rainfall is the only input of autogenic recharge in the karst unit, but the value of evapotranspiration should be discounted from the rainfall measured. This value is defined as effective rain. Jemcov and Petric (2009) assumed the effective infiltration as the input function instead of precipitation in order to take into consideration various processes in air, vegetation and soil and their influence on recharge, and to separate them from the processes within the karst aquifer system. In order to adequate the measured raw rainfall to the recharge area characteristics, the evapotranspiration has been taken account to get the effective rainfall. The rainfall (R) is used to soak up the overstorey, understorey and litter/moss zone (S) until a threshold (given for each type of soil) and for the evapotranspiration process (ETP). The remaining water is the percolated flow inward the unsaturated zone, i.e., the Effective Rainfall Re. The relation is given by the expression: Re=R−(S+ETP)(3.3) Evapotranspiration (evaporation + transpiration) is a term used to describe all the processes by which liquid water becomes vapor, from the surface of water, snow field or land (evaporation) or through plant stomata (transpiration) to the atmosphere. Therefore, the calculation of the actual evapotranspiration losses are necessary in the estimation of recharge due to precipitation. Direct measurement of evapotranspiration losses is difficult, because of that the interdependence between different components of the soilvegetation-atmosphere system hinders the proper estimation of the evapotranspiration losses (El-Baroudy et al., 2010). There are different methods to calculate the evapotranspiration (Brutsaert, 1982). Most measurements are indirect and spatially bounded. The choice of the method for the calculation of potential evapotranspiration in this thesis was governed by considering sufficient degree of reliability and meteorological information, required by the formula. To estimate the evapotranspiration in the karst analyzed for this thesis, the expression from Eagleman (1967) has been selected: ETP =f.C.W.Emaxp100 −Hr(3.4) 81 3.2 Karst background KARST AND GLACIAL HYDROLOGY where, Emax =6.1e(17.1T/(234.2+T)); f is a correction coefficient, its value is 0.0329 for the daily time step and 0.0329/24 for hourly time step; W depends on speed wind (for a calm wind is 0.8); Hris the relative humidity; and C depends on air temperature T (Table 3.1): Table 3.1: Coefficient C in the Eagleman’s expression Coefficient Condition 0.63 T <= 0◦ 0.63 + 0.024T 0 ◦< T<=21◦ 1.13 T>21 3.2.4.2 Output: Discharge A distinguishing feature of karst aquifers is that most of the groundwater is discharged through a small number of large springs. In some karst aquifers, the entire discharge is through a single spring. One or few of these springs are the base-flow discharge and are called underflow springs, while other springs flow only during periods of high discharge and are called overflow springs (or trop plein in French) (Worthington, 2003). A lot of small tributaries converge inside the aquifer, to form a flow path that drains toward the surface through a spring. So, the discharge from karst springs is a composite of all water moving through the aquifer. The spring (or springs), therefore, is an optimum location for measuring some parameters that are the outputs of the karst system. The common way to study the discharge of one spring as a function of time is the hydrograph of a single event, with a falling limb named recession curve that occurs after rainfall has ended. The integral of the curve of a hydrograph is the flow-mass curve that represents a hydrologic quantity of the aquifer. The shape of the recession limb of the hydrographs of karst springs was studied from the end of the XIX century. Maillet (1905) pioneered the work using a first-order LTI model approach for the drain of a single reservoir. He defined the recession coefficient alpha in the Eq. 3.5, where the discharge is proportional at any time to the water volume stored. Q(t) = Q0e−α(t−t0)(3.5) where Q0is the spring discharge at the beginning of the recession (t=t0), and αis the recession coefficient with dimension T−1. Subsequent authors introduced the technique of composite hydrograph recessions. More coefficients, sometimes empirical and without physical meaning, were developed in order to cope with the real hydrographs. Ford and Williams (1989a) present a summary on that issue. Many authors, e.g. Trilla and Pascual 82 KARST AND GLACIAL HYDROLOGY 3.2 Karst background (1974), Milanovic (1976), Torbarov (1976) or Milanovic (1981) work under the hypothesis of the existence of conduits of different diameter draining at different rates. However, this approach was not satisfactory because in karst hydrogeology, the application of numerical methods requires a specially tailored methodology (Palmer, 2002). 3.2.5 Characterization of karst aquifers The first step in any model study is to set up a schematic representation (conceptual model) of a real system. The karst modelling approach chosen for this work is the global model, which can be characterized as a posteriori modeling or black-box modelling, following the methodology of system identification. Models can be classified into these three basic models (white, gray or black). Nevertheless, differential equations, aquifer geometry, and a set of flow parameters under certain boundaries and initial conditions, are applied in the basic models to build a model adapted to a specific karst type (Wheater et al., 1993). Black box models, based in the convolution integral between the rainfall recharge and spring discharge are the simplest models for the study of karst hydrology. Black-box models are fed by the observations and with the aim of characterizing the catchment response from observed data, which determine the model structure. The mathematical relations between the input and output time series are derived without reference to physical laws. This makes them especially adequate in karst aquifers where the hydrogeological observation is limited by the great complexity and discontinuity of the medium (Padilla and Pulido, 1995). Black box model approach is often performed using time series analysis to study the correlation and spectral analysis of the karst system (Larocque et al., 1998). The input time series is the effective rainfall over the catchment, and the output time series is total discharge through springs. 3.2.6 Survey of techniques applied to karst Below, there is a review of the main techniques of system identification and analysis, revised in Chapter Chapter 2, than have been applied in karst systems. 83 3.2 Karst background KARST AND GLACIAL HYDROLOGY 3.2.6.1 Correlation analysis The strong and fast relationship between rainfall and the rising of discharge at typical karst springs was observed by local populations since ancient times. Analytical models of groundwater processes had been developed by the end of XIX century, followed by numerical methods based partially on physical conception. Given the heterogeneous nature of this process both in its spatial structure and in its temporal behavior, Mangin (1974), in his doctoral work, concluded that it is difficult to obtain a physical model that reproduces the expected spring hydrographs in response to a specified precipitation. Mangin suggested a longer time for the autocorrelogram of time-series discharge data and a value below 0.2 indicates greater aquifer inertia (system memory). A longer time delay between zero and the peak of a cross-correlation between precipitation and discharge time-series data indicates greater transit time of a signal within an aquifer. He works at the hydrographs of three classical springs located in Pyrenees: Aliou, Baget and Fontestorbes (Mangin, 1981b) Similar methods were applied to other karst springs mostly around the Mediterranean Sea and employed to characterize different types of karst (Freixes et al., 1996; Jiménez et al., 2001; Jiménez et al., 2002; Amraoui et al., 2004; Pérez et al., 2004; Valdes et al., 2006; Andreo et al., 2006). However, these techniques are well known in the system engineering, and they should be used very carefully following the engineering formalization of a system, with taking the output as the discharge and the input as the rainfall. 3.2.6.2 Cross-correlation Following the basic work of Box and Jenkins (1976), several authors introduced the cross correlation and spectral analysis relating input and output in hydrology field. Examples of those papers are: (Mangin, 1981a,9; Mangin and Pulido, 1983; Benavente et al., 1985; Cruz et al., 1987; Pulido et al., 1987; Mangin and Pulido-Bosch, 1991; Morales and Antiguedad, 1992; Rodríguez et al., 1994). A very complete study of correlation and cross-spectral analysis to characterize the transformations between the input function (precipitation) and the output function (discharge) was developed by Padilla and Pulido (1995). The parameters that can be deduced are the response time, the distinction between quickflow, intermediate flow and baseflow, and the mean delay. The method offers quantifiable and objective criteria for differentiation and comparisons of karst aquifers. The study of correlogram and amplitude cross functions turn out obtain the duration of the impulse response of the karst aquifer. The gain function is used to differentiate quickflow, intermediate flow and base flow. The delay in the responses can be deduced from the phase and coherence functions. 84 KARST AND GLACIAL HYDROLOGY 3.3 Glacier background Fig. 3.3: Earth’s annual and global mean energy balance from Kiehl and Trenberth (1997). 3.3.2 Glacier classification A traditional classification of glaciers according to size is: •Mountain or alpine glaciers are the smallest type of glacier. These glaciers can range in size from a small cirque to a larger development filling a mountain valley. •Piedmont glacier is a huge mass of snow formed by coalescence of many glaciers in a plain flat region at the foot of the mountain massif. Piedmont glaciers are between several thousand to several tens of thousands of square kilometers in size. •Continental glacier is the largest, with surface coverage in the order of 5 million square kilometers. Antarctica is a good example of a continental glacier. Another classification, usually discussed among geologists, is the glacier taxonomy according to the heat state and thermal regime of the ice layer. Alhamann (1933) considered the thermo-physical character of ice masses as a basis for differentiating glaciers. The glacier is tempered when its ice is at the temperature of fusion. Consequently, and based on his polar expeditions on 1930, Alhamann established two groups of glaciers: polar and temperate. About the same time, Lagally (1932) presented similar propose, sub-dividing glaciers into corresponding thermodynamic categories: kalt (cold) and warmen (warm). The temperature of a polar glacier (or cold), was perennially sub-freezing throughout, except for a shallow surface which might be warmed for a few centimeters each year by seasonal atmospheric variations. On the contrary, in a temperate glacier (or warm), the temperature is at the melting point except in the uppermost layers. 91 3.3 Glacier background KARST AND GLACIAL HYDROLOGY Since this thermodynamic connotation, glaciers of the polar type may exist at relatively low latitudes if their elevations are sufficiently great. At the same way, temperate glaciers may be found even inside the Antarctic/Arctic Circle at elevations low enough and sheltered from extreme chill weather. So, regardless of geographical location, the internal temperature average of the glacier can be used to identify it. Though, these implications are based on changing – and usually difficult to measure – thermo-physical characteristic. Alhamann (1933) faced this problem by introducing the sub-arctic type –later he called it sub-polar–, for glaciers where the seasonal warmth penetrates at a depth significantly greater than the one experienced in summer on polar glaciers and at lesser depth than the temperated type. By the same year, Lagally also recognized this intermediate class of glacier and called it transitional. To avoid some confusion from alternate application of these different terms, Miller (1976) suggested including both transitional categories with other two additional types: subtemperate and poly-thermal. To elucidate the characteristics of each of these five categories, Miller suggested thermal parameters as temperature limits, which are summarized in the table 3.2. Table 3.2: Classification of glaciers with the main characteristics, according to Miller (1976) . TYPE EXPLANATION THERMAL BOUNDS Polar Even in summer have negative temperature down a great depth and have no melting. All ice is entirely cold. -10 to -70 ◦C Subpolar Thickness negative temperature dominate but in summer melting is possible -2 to -10◦C Subtemperated to typify ice sheets during changes from fully polar to fully temperate 0 to -2◦C Temperated The whole thickness, except the uppermost layers, is at the melting temperature of the ice Poly-thermal The temperature range comprises at least two of the foregoing temperature zones Miller (1976) illustrates with examples the five-type classification of glaciers in the following areas: Antarctic and Greenland ice sheets (polar and poly-thermal), the Nepal Himalaya, Svalbard (polar to sub-polar), Lapland (sub-polar), sub-Arctic Norway (subtemperate), the Alps (poly thermal to temperate), the Canadian Rockies (sub-temperate), and the Juneau Icefield, Alaska (sub-temperate to temperate). Eraso and Domínguez (2006) claimed that glacier Collins, in the insular Antarctica and target of this study, is a temperated glacier. 92 KARST AND GLACIAL HYDROLOGY 3.3 Glacier background 3.3.3 Glacier hydrology From a hydrological perspective, glaciers represent important water resources, contributing significantly to streamflow. Glaciers exert a considerable influence on catchment hydrology, by temporarily storing water as snow and ice on many timescales (Jansson et al., 2003). Understanding water movement through a glacier is fundamental to understand the glacier dynamics, induced floods, and the prediction of runoff throughout drainage basins. Nevertheless, glacier hydrology is a set of complex processes, since the released water suffers several changes of state and comes from several origins. The physical structure of the snowpack is a porous matrix, and will be able to hold the liquid water between the snow grains and will increase snow density. If the snowpack voids, where the liquid water is retained, become fill up, some of the snowmelt will start to fall down and take part in some of the glacier flows. Additionally, some of the snowmelt can more deeply infiltrate into the ground to form the groundwater flows. A simplified hydrological process development of a glacier hydrology is depicted in the longitudinal cross section of Fig. 3.4. After rain falls or snowpack melt, water can run along the glacier surface (supraglacial osurface runoff), percolates through snow by crevices and moulins to form the endoglacier flow, or remains dammed. An amount of this water can reach the glacier bed to set a subglacier flow or be retained in an inner resevoir. Finally, part of liquid can move down-gradient till the water table, to form groundwater with the common features of a hydrology underground. A complete knowledge of a glacial hydrology comprises the total contribution of glacier surface runoff, englacial water storage and transport, subglacial drainage, and subsurface groundwater flow. The four components are represented by Flowers and Clarke (2002), as a two-dimensional model with vertically integrated layers that exchange water among them. This model assumes that incompressibility and water volume are conserved according to Eq. 3.6, after simplifying some interglacial flows: ∂hr ∂t+∂Qr ∂x=M+R−φr:e+φr:s−φr:a(3.6) where r superindex indicates the surface origin of water, M is the melting rate, Qr is the water discharge per unit width (or flux) due to superficial runoff, R is the rate of precipitation falling as rain, φr:eis a source/sink term for water exchange with the englacial layers, φr:sis the rate to subglacial, φr:arepresents the exchange with the groundwater system. 93 3.3 Glacier background KARST AND GLACIAL HYDROLOGY Fig. 3.4: The hydrology development and melting process in a glacier, modified from Price et al. (1979) Some authors (Derikx, 1969) simulate the discharge from a glacier by considering the analogy with a ground-water system. Hock and Janssen (2006) observe two spheres of activities in studying glaciers. Though the boundary is diffuse, from the front of the glacier upward could be the general space for glaciologists, whilst from the front of glacier downward is rather an area for hydrologists. Nevertheless, the common input for both is the weather phenomenon, and under the glacier mass, the groundwater movement is similar to a karst aquifer. Fig. 3.4 depicts unsaturated and saturated zones collecting water from glacier melting and rain percolation, as a karst aquifer receives water percolated from the rainfall. Both systems retain water and deliver it due to different circumstances, being the air temperature the main cause of the drainage of the glacier. This similarity in cryokarst term was assumed by Eraso and Domínguez (2005). Assuming these parallel concepts, some data-driven models widely applied on karst aquifers, could be tested on glacial hydrology, with the caveat of taking in mind the glacier system singularities. 3.3.4 Modeling approaches The aim of models should be predicting the glacier discharge as the integral response of a glacier basin to the variable weather conditions and to show the advance or retreat of the glacier front. It is unworkable correctly calculate the impact of all the geophysical 94 KARST AND GLACIAL HYDROLOGY 3.3 Glacier background variables and meteorological factors in time and space. Therefore, a careful study of the most significant input variables of the system, could simplify the model with acceptable results. Basic model types are classified according two main categories (Fig. 3.5): Fig. 3.5: Impact weather variables on glacier discharge, categorized by melting models and runoff models. a) Melting models. There is a complete hierarchy of melt models relating ablation to meteorological conditions, varying greatly in complexity and scope. Classically in glaciology two main different approaches can be distinguish in snowmelt simulation. One is well grounded in the detailed evaluation of the surface energy fluxes (energy-balance models), based on the solution of the equation of energy balance of the snow fallen and melted. The second class of models considers the weather variables as indicators of physical processes. Concretely, the air temperature is used as unique index of melt energy (temperature-index models and degree-day method). The degree-day method has been applied in many cases for more than a century. However, these are often not practical due to large data requirements and uncertainties about spatial variability. The temperatureindex method, most commonly used, generally provides better performance by the low data requirements and simplicity. b) Runoff models. Among runoff models, some are rather linked to physical structures by setting equations derived from the law of mass conservation, and expressed as a balance between the internal distribution of water and external sources, e.g. Flowers and Clarke (2002). Others are based on weather variables –air temperature, precipitation or 95 3.3 Glacier background KARST AND GLACIAL HYDROLOGY radiation– as an input and the discharge of glacier-fed streams as output. The model is developed from actual daily observations of discharge and simultaneous weather conditions, taking the dynamics of glacier as a black-box system. These models provide practical solutions to hydrological problems, because of using an analytical method under a low computational cost, and requiring simple information, though long and dense time series are preferable. An example is a conceptual model identified with a transfer function or other analytic model structures, established through observations of input and output relations and constrained to a set of boundary conditions. 3.3.4.1 Input variables The main posed problem was to set the weight of each meteorological events which participate in the melting process. Some classic discussions have emerged to resolve which of the two phenomena, air temperature or solar radiation, have a more direct influence on the melting. Walcher (1773) was one of the first to propose that glacier fluctuations are caused by variations in climatic conditions. Later, Mousson (1854) mentions the effects of solar radiation, the air contact, and rainfall as the three ways in which ablation energy is delivered to the glacier surface. Finsterwalder and Schunk (1887) already assumed a relation between ablation and temperature for their study of the variation of the Suldenferner glacier in the Eastern Alps. A little later, Hess (1904) clearly stated that radiation is the most important cause of ablation. Numerous studies, gathered by Paterson (1994), have been carried out from 1930 to measure the components of the glacier surface energy balance. In most of them the net radiation constitutes more than half of the total energy supply. Nevertheless, Angström (1933) stressed the importance of several agents as temperature, radiation and wind as agents for melting. Subsequent papers have compared mass-balance with the air temperature data taken from weather stations close to glacier, although in some cases also consider the precipitation as a part of the input. Among them, Wallen (1949) thinks that there is a moderate correlation between mass balance and air temperature vary. Alhamann (1953) highlights the length of the period when the temperature is above the melting point, in addition other factors as the annual precipitation, the amount of incoming and outgoing radiation influenced by the degree of cloudiness, wind velocity and humidity. Hoinkes (1955) wrote: In recent years many authors, on the basis of careful studies, have come to the conclusion that summer temperature is to be regarded as the most important factor influencing the behaviour of glaciers. Also Hoinkes (1955) explicitly 96 KARST AND GLACIAL HYDROLOGY 3.3 Glacier background rejects that the heat exchange from air to ice during a hot summer is sufficient to account for the greater ablation. Many recent studies have revealed a high correlation between melt and air temperature. Braithwaite and Olesen (1988) confirm that net radiation is the main energy source for melt, but the air temperature is usually better correlated with melt production and run-off than net radiation. The air temperature is used as a proxy for melt energy input to the glacier, because it has a strong relation with solar radiation (Lafrenière et al., 2003; Jansson et al., 2003). Braithwaite and Olesen (1988) found a correlation coefficient of 0.96 between annual ice ablation and positive air temperature sums in some Alpine glaciers. Then the air temperature may act as the sole index of melt energy in spite of the predominance of net radiation. It is attributed to the high correlation of temperature with several energy balance components (Ambach, 1988; Kuhn, 1993). 3.3.4.2 Glacier discharge Glacier melt provides important contributions to surrounding areas through the streamflows (or glacial streams) that surge at the front of the glacier. Reversely, streamflows are a reflection of input variations of the glacier and its transformation in storage and transmission along the glacier course. Glacier melting oscillations is reflected on water discharge. The runoff has a marked diurnal peak due to the diurnal radiation cycle, with a typical delay of several hours between peak melting and peak at discharge. The annual cycle is divided between a cold and inert period, with snow accumulation, and a warmer, active period. Between the cold and warm periods, and vice-versa, there are transition seasons with complex features. The beginning is dominated by the melting of the accumulated snow, with poor water flow and refreezing processes. With the advance of the season, the melting area increases and the ice starts to melt when it is uncovered. The process continues with the formation and/or reactivation of a number of preferential vertical and horizontal flow paths on, in and under the ice, through crevasses and horizontal conduits. The melting water, routed from supraglacial, reaches the englacial and sub-glacial flow, and finally becomes the proglacial discharge (Fig 3.6). All the hydrological gauges and assessments of the water discharge, located after the front of glacier, can be adequate when the direct measurements and observations of glacier mass balance are difficult or perhaps impracticable. The streamflow hydrograph is influenced both by the snow mass availability and by seasonally variations in the weather factors. So, the analysis of trends in magnitude and timing in a streamflow 97 3.3 Glacier background KARST AND GLACIAL HYDROLOGY provides sufficient and confident indicators for trends and fluctuations of the whole glacier discharge, and consequently 0of climate changes. 3.3.4.3 Energy balance model (EBM) Knowledge of a glacier’s energy budget is important to estimate the contribution of climatological factors to glaciers. The problem should be posed on how the free atmosphere interacts with the glacier surface. Glacier-atmosphere energy balance can be explained from the exergy concept of a system. Exergy is defined as the maximum work which can be produced by a system, a flow of matter or energy, to reach the thermodynamic equilibrium with a specified reference environment. Exergy inflows, from the sun to the earth, must be balanced by exergy outflows (in the form of IR heat radiation). The exergy terrestrial inventory is varying, depending on biospheric processes and antropogenic activities. The difference between solar exergy inflows and Earth’s exergy outflows is the energy available to develop work, i.e. to carry out all atmospheric, hydrospheric, geospheric and biospheric processes. Following the expression for exergy given by (Valero and Valero, 2010), the change of exergy for the Earth is (3.7): ∆B=Z2 11−T0 T1∂Q−[W−p0(V2−V1)]−T0σ;[in joules](3.7) where T0and p0denotes the temperature and pressure at ambient conditions and Tjis the surface temperature where the heat transfer takes place. The integral term represents the exergy transfer accompanying heat, the term [W−p0(V2−V1)] is the exergy transfer accompanying work, and the term T0σaccounts for the destruction of exergy or irreversibility. The hydrosphere extracts exergy from the solar inflow (in the form of heat or work) in order to accumulate it in somehow, or consume it in physical processes, as the ones come about in glaciers. So, snow and ice play an important role, interacting with the atmosphere over a range of spatial and temporal scales in a complex feedback mechanisms. Once calculated the whole exergy of a coupled system formed by climate and glacier subsystems, the exergy component due to transfer accompanying heat is the responsible of governing the melting process and consequently the discharge (Fig. 3.6). An imbalance in the thermal exergy is an available energy budget that will affect the temperature of the glacier surface. If the glacier surface is below but near the melting point (0 ◦C), an additional contribution of heat transferred towards the glacier surface, will yield the melt of ice or snow. The sources energy that cause snow melt are the sum of 98 KARST AND GLACIAL HYDROLOGY 3.3 Glacier background Fig. 3.6: Diagram illustrating that the exergy component due to transfer accompanying heat is the responsible of governing the melting process and consequently the discharge shortwave Qsn net radiation, long-wave net radiation ln, convection from the air Qh, vapor condensation Qe, conduction from the ground Qg, and energy contained in rainfall Qp. The available energy budget Qmis expressed by the equation (3.8), given by Price et al. (1979): Qm=Qsn +Qln +Qh+Qe+Qg+Qp−∆Qi(3.8) These contributions are usually measured as energy per time per unit area of snow. ∆Qiis the rate of change in the internal energy stored in the snow per unit area of snowpack. In the warm period, the net flux of heat (∆Qi) goes inwardly the snow, while during periods of cooling, the net flux (∆Qi) emerges outwardly from the snowpack. If the available energy budget is sufficient, the melting process can produce liquid water. The Fig. 3.4 is a graphical synopsis of the melting process, where the energy flux is representing by red lines. Many authors have studied the energy and mass balance of glaciers. Pioneering work concerning the details of energy exchange between the atmosphere and a glacier surface was performed in the Nordic countries. Alhamann (1948) derived the first empirical formula for the computation of ablation from known values of incident radiation, air 99 3.3 Glacier background KARST AND GLACIAL HYDROLOGY temperature and wind velocity. Sverdrup (1935) computed a complete energy balance. In subsequent decades, gradient flux techniques to ice and snow were studied (Ambach, 1963; Munro, 1990; Hay and Fitzharris, 1988). Others consider the mass balance depends strongly on altitude, (Braithwaite and Olesen, 1988; Munro, 1990; Greuell et al., 1997). Others have modeled mass balance gradients and sensitivities (Ambach, 1963; Oerlemans, 1992; Van de Wal and Oerlemans, 1994; Jóhannesson, 1997). Only a few authors have studied the surface of a glacier in a three-dimensional way (Arnold et al., 1996; Oerlemans et al., 1999). 3.3.4.4 Temperature-index model (TIM) This method computes snowmelt only by using air temperature as an index to melting process. A thermal index is a statistical estimator used instead of physical variables. This statistical method removes many of the variables from the equations, because is often difficult to assess. Moreover, air temperature is assumed as a predominant variable in the energy budget equations, because it is connected with many of the energy exchanges involved in snowmelt. The basic equation for the temperature index solution is defined in (3.9), although Gray and Prowse (1993) provide some derived expressions. Ms=Cm(Ta−Tb)(3.9) where Ms= snowmelt, in. per period, Cm= melt-rate coefficient that is often variable, in./(degree/period), Ta= air temperature in ◦T, and Tb= base temperature in ◦F. Degree-day model is a form of temperature index models that are based on an assumed relationship between ablation and air temperature usually expressed in the form of positive temperature sums. Hock (2003) relates the amount of ice melt, M (mm) in The Eq. 3.10, during a period of n time intervals, dt (day), to the sum of positive air temperatures of each time interval. n ∑ i=1 M=DDF n ∑ i=1 T+∆t(3.10) The factor of proportionality is the degree-day factor (DDF), and is expressed in mm d−1 ◦C−1. Commonly, a daily time interval is used for temperature integration, although any other time interval, such as hourly or monthly can also be used for determining degreeday factors. 100 Chapter 4 Analysis and Identification of Fuenmayor aquifer ... a country like ours, where, because of its elevated heights average, its rivers have to pour waters tumultuously, in a country so unfortunate as ours where the cries of pain by floods drown the anguish caused by drought... Lucas Mallada (1841-1921) 4.1 Introduction Fuenmayor spring has been monitored continuously for identification purposes to study the behavior of a karst groundwater system. Under a linear time invariant hypothesis, the application of the simple correlation of spectral analysis and parametric identification of the transfer function generated some interesting results in the monitored spring. These tools have historically been successful in studying a large number of karst springs and continue to be practical approximations in initial attempts to obtain a draft model. Because of the nonlinear and nonstationary nature of karst, more effective systemic techniques are required to cover certain aspects of analysis that the linear system cannot reveal adequately. This chapter presents interesting results using Fuenmayor spring data, collected over several years, to show the ability of wavelet techniques in the analysis and identification of a karst spring system. The essential content of this chapter has been published by Chinarro et al. (2011) and Cuchí et al. (2013). The objective of this chapter This chapter aims the system identification and analysis the time series of effective rainfall and discharge in Fuenmayor spring, trying diverse structures in order to define linear and non linear models that characterize karst hydrology. 107 4.1 Introduction FUENMAYOR AQUIFER Outline of this chapter This chapter is organized as follows: Sec. 4.2 is the geographical framework of the area under study in Natural de la Sierra y Cañones de Guara (Natural Park of Guara Mountains and Canyons) located in the central and highest part of the Pyrenean External Sierras, where Fuenmayor spring is situated. An overview on Fuenmayor spring is in Sec. (4.3), with a description of karst system and how the water appears at the contact between the limestone of the Guara formation and bounding of recharge area around hills surrounding Ciano polje (Sec. 4.3.1). A hydrological overview is given in Sec. 4.3.2, where the behavior of the aquifer is known by the evolution of its spring discharge and by a strong relationship between rainfall and discharge. Also, a literature review on this spring are depicted. Sec. 4.3.3 is a list of infrastructure and facilities that have supported the monitoring and exploitation of the aquifer hydrology data. Section (4.4) is an analysis of signals or time series study. A paragraph highlights the importance of choosing hourly sampling period (4.4.1) in order to avoid the loss of information, and issues about preparation data (4.4.2). In order to find causality between the precipitation phenomena and the aquifer discharge, Sec.4.4.3 presents simple correlation and spectral analysis with results about autocorrelation and the spectral density of the discharge of the Fuenmayor spring. In addition, Sec. 4.4.4 uses wavelet power spectrum to reveal cyclic behavior patterns in precipitation and discharge signals. The analysis in this section assumes the nonlinearity of this karst spring system to discuss the results from scalogram observations, proposes a method for rainfall-discharge model calibration and and estimate the representation of cross power spectrum and coherence. The key problem in the system identification is to find a model structure with sufficient flexibility to suit the karst dynamics. System identification procedures start with nonparametric identification, showing the hydrogram and recession limb (Sec. 4.5.1.1), going on by creating model based on cross-correlation and impulse response. The Parametric Identification of Fuenmayor section deals with the implementation of linear methods (4.4.3) for the study of the behavior of this spring, to obtain the first transfer function in this aquifer, proposed under the hypothesis of a linear and time invariant response. Using previous results from linear identification and spectral analysis, the final stage in this chapter consists of tuning a nonlinear model from HammersteinWiener structure. Finally, in the Conclusion section, there are important assertions related to linear and non-linear identification of karst springs, and how to interpret the coherence results. 108 FUENMAYOR AQUIFER 4.2 Geographical framework 4.2 Geographical framework PARQUE Natural de la Sierra y Cañones de Guara (Natural Park of Guara Mountains and Canyons) is located in the central and highest part of the Pyrenean External Sierras, forming a large limestone barrier that constitutes the southern limit of the Pyrenees of Huesca, Spain. The site, declared Natural Park by Courts of Aragon in 1990, has a surface of 47,450 hectares, in addition to 33,775 hectares of Peripheral Protection Zone, involves 15 municipalities, and comprises the Sierra de Guara, Gabardiella, Arangol, Balcés and Sevil. The Tozal de Guara, with 2.077 m of altitude, is the highest point of the park. Fig. 4.1: Huesca province containing Parque Natural de la Sierra y Cañones de Guara. The green triangle is the location of Fuenmayor spring. Map from Cuchí and Setrini (2004). Guara mountains represents the southern part of the mantle of Gavarnie. The most important geological materials are the fissured and karstified Eocene limestone, reaching thicknesses close to one thousand meters. There are, nevertheless, several impervious units delimiting diverse hydrogeological units (Puyal et al., 1998). The thickness of Guara 109 4.3 Fuenmayor spring FUENMAYOR AQUIFER Formation limestone increases from west to east and is intensely fractured, which favors its dissolution. The alpine orogenic creates a series of anticlinal and synclinal structures that have provided the establishment of major structural forms and a dense network of limestone fractures. (Cuchí, 1998). During the last millions of years, an intense process of karstifications has been performed in Guara Mountains. The most important karst features are the doline field of Cupiarlo, and the poljes of Ciano and Vallimona (Rodriguez-Vidal, 1986). The hydrology is directly related to karst features. The area is drained by four major rivers: Flumen, Guatizalema, Alcanadre and Vero, deeply incesed by stunning canyons (Fig. 4.1). The underground of the Guara Mountains shows sinkholes and caves produced by groundwater flow under karst conditions. In the area, several karst aquifers have been identified by their springs, for example, Trinity of Rasal, Petrolanga, Solencios of Bastarás, Santa Cilia, Morrano and Pedruel, and Fuenmayor spring. The last one is the focus of this study. The area has a typical mediterranean weather with rainfall periods in Spring and Fall. The average rainfall, from 1955 until 1984, was 793 mm. The interannual precipitation shows the typical irregularity of the Mediterranean climate where droughts are common. The average temperature is near 12 degrees Celsius. 4.3 Fuenmayor spring The karst spring of Fuenmayor, located at 714 m (a.s.l.), emerges near San Julian de Banzo (Fig. 4.2.b), which is a small village with about 50 inhabitants, located at the north of Barluenga and Chibluco, in the province of Huesca (Spain). The ravine of San Julian comes from the spring and flows into the Flumen river; it is the axis of the Fuenmayor basin and natural separation of the two districts of the town, Suso and Yuso (fig. 4.3). Fuenmayor spring appears at the foothill of the External sierras of the south central Pyrenees (Millán, 1996). Located about 20 km from the city of Huesca, spring has been used for drinking water supply in the city since 1880, the year in which a cholera epidemic propagated in the area near the fluvial resources which forced the seeking of other alternatives (Cuchí et al., 2002). 110 FUENMAYOR AQUIFER 4.3 Fuenmayor spring (a) Polje of Ciano, North of San Julian de Banzo. It is the northern area of Fuenmayor aquifer recharge A canyon runs dry from the small Ciano polje until the spring (b) The headwater basin of Fuenmayor spring (a) and an overview of San Julián de Banzo, with Salto de Roldan canyon on the background (b) Fig. 4.2: An overview of San Julián de Banzo town and basin of Fuenmayor spring. 4.3.1 Geological description San Julian de Banzo village, located at the startup of glacis Somontano, is dominated by the escarpment of Pre-Pyrenees conglomerates formed by Salto del Roldan canyon(1121 m) and San Martin de la Baldosera (1431 m), and limestone of Serreta de Vallés (1128 m) and Matapaños (1532 m). The southern boundary is formed by sandstone and clay of Sariñena formation (Miocene), covered on terraces and glacis. Structurally the area is part of the inflexion in Southern Pyrenean frontal overthrusting, between San Julian de Banzo and La Almunia del Romeral. It consists of a series of folds with an east-west orientation (Cuchí et al., 2006). Eocene and Upper Cretaceous limestone host a modest karst aquifer, drained by Fuenmayor, located between the back overthrusts of San Julián and Cuello Bail (Fig. 4.4). Fuenmayor is a spring of dammed type, without trop pleins. The water appears at the contact between the limestone of the Guara formation (Middle Eocene) and the sandstone of the Sariñena formation (Miocene). The spring is located in the middle reach of the Molon creek. Its small canyon runs dry from the small Ciano polje (Fig. 4.2.a) until the spring. Local population believed that recharge area of the spring is located in the hills around Ciano polje. Trilla and Pascual (1974) estimated a recharge area under 15 km2. This value was agreed by Pinilla (1996). The recharge area was also calculated by Cuchí et al. (2002), following a structural and hydrological criteria, to be near 10 km2(Fig. 4.5). Later, Oliván (2013) points out an area of 9,70 km2. An extra allogenic recharge may arrive from the San Martin conglomerates through the homonomous creek (Cuchí and Villarroel, 2002). In the associated recharge area, the Keuper clays acted as the impervious base (Pinilla, 1996). 111 4.3 Fuenmayor spring FUENMAYOR AQUIFER Fig. 4.3: Aerial view of the source area of Fuenmayor spring, marked by a red circle. In both side of the stream, in the left corner of the picture, are the districts of San Julian de Banzo (Suso to the left, and Yuso to the right). Courtesy of Google-Map. Near Fuenmayor there is Dos Caños spring situated about 500 meters southeast of the Fuenmayor spring (Fig. 4.5). Dos Caños drains a small aquifer with a recharge area does not exceed the square kilometer. Initially Dos Caños was considered independent of Fuenmayor in their hydrology and hydrochemistry data. Water from the Dos Caños is warmer and saltier than the water from Fuenmayor. Both aquifers are limited by the supposed impervious clays and evaporites of the Keuper facies (Pinilla, 1996). However, a pumping test carried out in the only well tapping Dos Caños in 2005 also affected Fuenmayor spring, indicating a hydrological connection, suggested by Cuchí and Villarroel (2002). These authors proposed two hypotheses about interconnection between both aquifers. First one, there are approach points to a shallow underground connection of both aquifers. Another hypothesis establishes a much deeper connection that separates two bodies of different salinity by an interface (Cuchí and Villarroel, 2002). According to Cuchí et al. (2006), the recharge area shows clear exokarst features. Lapiaz structures, mostly hohlkarren type (Sweeting, 1973), are formed by disolution of the surface of the limestone. No karst caves have been discovered in the area. A large part of the recharge area is covered with lithosol (like a barren area). The rest shows thin 112 FUENMAYOR AQUIFER 4.3 Fuenmayor spring Fig. 4.4: Geological scheme of the study area, according to Millán (1996) soil. Maximum humidity storage of the soil has been estimated to be 30 mm (Cuchí et al., 2006). 4.3.2 Hydrological overview Fuenmayor spring already was referred by Mallada (1878). The behavior of the aquifer is known by the evolution of its spring discharge. From early times, it is known that Fuenmayor presents a strong relationship between rainfall and discharge, with short floods after raining and long low discharge periods during the summers. For this reason, the karst nature of the aquifer was recognized by Trilla and Pascual (1974). Consequently, several actions have been performed in order to regulate the spring, as the building of a water mine suggested in 1955. Later, other actions, as those summarized by Cuchí et al. (2006), were made at Fuenmayor and surroundings, in order to increase the water supply to the city of Huesca. For this purpose, a unsuccessful borehole was drilled at the spring in 1984. Because of this failure, at 1990, two manual water gauges were built by the IGME (Geological Survey of Spain) and CHE (Water Authority of Ebro Watershed) in order to obtain more information on the spring behavior. One was built at the Huesca aqueduct and the other at the spillway of spring. 113 4.3 Fuenmayor spring FUENMAYOR AQUIFER Fig. 4.5: Schema of two structural hydrogeological units: Fuenmayor (center), with its hypothetical basin limits, and Dos Caños (south) (Cuchí and Villarroel, 2002) on geological section. At 2000, the GTE, Grupo de Tecnologías en Entornos hostiles (Group of Technologies in Hostile Environments), of the Zaragoza University designed and implemented an automatic monitoring station to measure continuously the following parameters: discharge, rainfall, electrical conductivity, water and air temperature. The device was working uninterruptedly until the summer of 2005, when a pumping test from one, until that moment, unused well, dried the spring. A literature review of work carried out on Fuenmayor spring is collected in Table 4.1. Table 4.1: Literature review of studies about Fuenmayor spring. DATE FACT TO BE MENTIONED CITATION From ancient times From ancient times, the water from Fuenmayor is used for irrigation and water supply to the local population. Cuchí et al. (2002) 1878 Description of this spring in the publication about the geology of Huesca province. Mallada (1878) (Continued on the next page) 114 FUENMAYOR AQUIFER 4.3 Fuenmayor spring Table 4.1: (Continuation) DATE FACT TO BE MENTIONED CITATION 19th century Construction of the water conduct to the city of Huesca as an alternative to the use of springs of Ibón and Ángel that capture waters of Huesca watershed. Cuchí et al. (2002) 1955 The first one who estimated the water balance of Fuenmayor spring and proposed the construction of a drainage gallery. Lasierra (1955) 1974 They explicitly recognize the karst nature of the aquifer, initiate systematic appraisals, correlate discharge with rain and analyze the flow depletion curve. Trilla and Pascual (1974) 1984 City of Huesca performed by percussion, a survey of 60 m. in the Gorga Mora, at the header of the spring but it does not have the expected success. Villarroel and Cuchí. (2004) 1990 Confederación Hidrográfica del Ebro installs two gauging weirs in the output of spring, in the Huesca pipe and the spillway of the spring. Initially measurements are controlled visually. Villarroel and Cuchí. (2004) 1992 A survey for research was performed at level 787, by direct rotation, in the vicinity of Fuenmayor, and 237 meters were tested (report 2911-7-0012). Carceller (2007) 1994 At about 800 m to the southeastern spring, a pilot drilling was made by percussion, reaching 300 meters deep. This survey captured an aquifer, that drains by the Dos Caños spring. (report 2911-7-0013) Carceller (2007) 1995 Publication on the hydrogeological potential of the thrust structures on the southern edge of the Sierra de Guara. Octavio et al. (1995) 2000 The research group of University of Zaragoza, GTE, installs a continuous monitoring station of output flow, precipitation, electrical conductivity and temperature of water and atmosphere. Villarroel and Cuchí. (2004) 2003 Water-rock interaction in the referred area. Monaj (2003) 2011 Application of cross wavelet power spectrum to quantify an indication of the similarity or strong relation of rainfalldischarge, independently of the energy amplitude of the signals, and to explain some hydrological cycle phenomena in the karst system. Chinarro (2009) 2011 About linear methods employing system engineering techniques for the analysis and identification of a hydrological system, considered as a black box model, which disregards information on the internal structure of the aquifer. Under a linear time invariant hypothesis, the application of the simple correlation of spectral analysis and parametric identification of the transfer function generated some interesting results in the monitored spring. Chinarro et al. (2011) 2011 Hydrogeochemical and isotopic characterization of karstcarbonate aquifer. Oliván et al. (2011) (Continued on the next page) 115 4.3 Fuenmayor spring FUENMAYOR AQUIFER Table 4.1: (Continuation) DATE FACT TO BE MENTIONED CITATION 2013 Identification techniques used to study the relationship between rainfall and the discharge of karst aquifers in order to examine some of the possibilities they offer and to address their limitations. The paper applies the reviewed techniques to a time series of the Fuenmayor karst spring, in the southern central Pyrenees. An acceptable linear response is shown. The quality of the different models obtained is evaluated with the Nash-Shutcliffe model efficiency coefficient. Cuchí et al. (2013) 2013 Delimitation, assessment of recharge, and aquifer hydrodynamics drained by karst spring Fuenmayor Oliván (2013) 4.3.3 Infrastructure and facilities The outflow of Fuenmayor karst system flows, from north to south, in an open stream in about 30 m to reach a diversion dam. It is constructed of concrete, with an embankment of 12 m long, a depth not exceeding 0.5 m and a covering area of 50 m2. Nearby, a small squared building of 10 m2provides sheltering for maintenance tools and monitoring equipment. In the basement, water flows are channeled through a gate to manually regulate the flow toward two possible outputs. The first one is to supply piped fresh water to the city of Huesca, with a Thompson weir to measure the water flow. The other output is a bypass to the natural bed of the spring outflows. Downstream, there is a second rectangular gauge to measure the flow from this channel and water from spillways of the dam. The hydrological monitoring station was installed in 2000 year, with the aim of measuring, on an hourly basis, discharge, accumulated rainfall, electrical conductivity of the water, and air and water temperature (Monaj, 2003). The main elements of this station are: •The output discharge is measured simultaneously at two outlets by weirs, with a pressure probe PDCR 130/D. They are pressure transducers with a range of measures from 70 mbar (millibar) to 135 mbar, with a linearity of +/- 0.1 % and a total error of +/- 1.5 % between -20 ◦C and +80 ◦C. It is an active device that requires a 12 v power supply. •A rain gauge, model 52202 Young, was attached on the roof of the building. It has an area of 200 cm2, captures with a resolution of 0.1 mm lever tips, and achieves the accuracy of 2 % in measurement up to 25 mm/hr and 3% up to 50 mm/hr. 116 FUENMAYOR AQUIFER 4.4 Analysis of signals analysis, the set of used scales can be chosen arbitrarily. By visual inspection of Fig. 4.10, some features are perceptible: •The spectral map shows a semiannual cycle following the typical pattern in the Mediterranean climate, which is observed by the continuous dark red in horizontal direction, in the scales range of 171-241. •Another highlight spot is located beyond the scale 341, closer to the scale 360, in all timeline values, which corresponds to the annual hydrological cycle. •The black dashed line represents the COI, to highlight the influence of the ends of the range of data. Since COI is a trust boundary, findings outside the cone would be under suspicion of unreliability. •In rainy seasons (6-13, 16-20, 23-25 month intervals), the rain reveals a richer spectrum with frequency components in the 200-2000 hour band (about 8-80 days) . In dry intervals, these components are absent. •The plot shows rainfall events as small isolated zones in the high frequency band that requires to zoom the spectral picture. 4.4.4.2 Wavelet power spectrum of discharge signal Analogously, at a glance on Fig. 4.10 down, the power spectrum of the output spring reveals also natural oscillations of discharge regime. The distribution of energy in the discharge signal and the concentration in certain areas can be observed by different color gradient levels. The hue gives a quantitative idea of the frequency contained in a signal at a certain time, so it represents the possible value range of frequencies contained in the discharge signal. Dark red peaks of the discharge emphasize visibly the annual and semiannual cycles. The wavelet window, acting as a filter, can be moved both along the time line and along the scales with a defined granulation. So, more details about a specific area of the spectral map can be obtained. The upper spectrum is zoomed, to observe, in high frequency bands, details on possible cycles of the discharge in short periods. Similarly, in case of an extended series of data, the filter can be shifted towards lower frequency bands (upper scales), to find likely rainy cycles and drought cycles along lengthy periods of time (Fig. 4.11). The longer the data series is, the less COI will influence the ends, and the more accuracy of correlation will be achieved, especially at low frequencies. 123 4.4 Analysis of signals FUENMAYOR AQUIFER Fig. 4.11: Wavelet spectral power of the rainfall (up) and the discharge (down) in Fuenmayor spring zoomed to display details in high frequency bands. It is interesting to know is that the wavelet transform can depict and place the assessment of the frequency component at the time it occurs. However, Fourier procedures can reveal that there are some certain frequencies in the signal, but not when they happen. By a simple visual examination of the Fig. 4.11, the wavelet spectral power exhibits a quasiperiodic behavior in the high frequency. The discharge power spectrum discloses the weight of each frequency carried by the signal. Between the dark blue areas that represent the common dry season of each year, there are a few spectral power peaks in the frequency range approximately between 30 and 43 days. This happens every year, although with slight differences of hue. These results are similar to those obtained by linear methods. The plot also shows erratic rainfall events forming small isolated zones in the very high frequency band. Obviously, dry intervals are characterized by the scarceness of spectral components. Because of the rainfall effects over discharge are attenuated in wet period, the karst groundwater system behaves as a loss-pass filter. Importantly, in the discharge 124 FUENMAYOR AQUIFER 4.5 System identification spectrum graph, the weaker color in certain areas, outside the COI bounds, indicates the area of unreliable values. Therefore, the resolution advantage of the wavelet transform permits detection of changes in nonstationary of discharge and rainfall, with exceptional localization both in time and frequency domains. 4.5 System identification The key problem in the system identification is to find a model structure, with sufficient flexibility to suit the karst dynamics. Once the signals observed in Fuenmayor spring have been analyzed, in the following paragraphs, the relationship between precipitation and discharge signals of Fuenmayor spring are studied following the different techniques reviewed in Chap. 2. 4.5.1 Nonparametric system identification This section applies the main nonparametric system identification methods, given in Sec. 2.7, to the effective rainfall and discharge os Fuenmayor spring. Fig. 4.12: cross-correlation between effective rainfall and discharge in Fuenmayor spring 125 4.5 System identification FUENMAYOR AQUIFER 4.5.1.1 Cross-correlation. Impulse response The cross-correlation coefficient between the effective discharge and discharge of the Fuenmayor spring is depicted in Fig. 4.12 following the procedures indicated in Sec. 2.7.1.1. An estimate of the unit hydrogram with a sharp pointed peak at 14 hours suggests a fast flow into a highly karstified system. The recession limb can be divided into two sections. The first one, from 14 hours to approximately ten days, has a fast decay. In the second, the decay is slower and more extended in time. This suggests that the behavior can be identified by a system of second order. Fig. 4.13: Kernel estimation of Fuenmayor spring based on cross-correlation: (a) without prewhitening and (b) with prewhitening. Using theory described in Sec.2.7.1.1, the cross-correlation coefficient to estimate the impulse response is represented in Fig. 4.13 a.1, and, applying the convolution sum, it is possible to simulate the system. Fig. 4.13 a.2 shows the result of this simulation for the effective rainfall in comparison with measured discharge. As can be seen, the discharge is greatly overestimated. This model based on the cross-correlation coefficient gives an efficiency factor of E = -314.14 for the Fuenmayor data applying Eq. 2.48, which gives a poor efficiency value, according discussed in Cuchí et al. (2013). However, as has been stated in section 3.2.4.1, the effective rainfall of the Fuenmayor karst system is far from being white noise. To overcome this problem, the whitening filter technique has been applied. Fig. 4.13 b.1 shows the impulsional response estimate obtained after the whitening. Fig. 4.13 b.2 shows that in this case the simulated discharge 126 FUENMAYOR AQUIFER 4.5 System identification Fig. 4.14: Linear coherence function of Fuenmayor spring fits much better than the previous one. The simulation model presents an efficiency value of E = 0.095. Finally, the estimated impulsional response is near zero after 100 days. Fig. 4.14 shows the linear coherence function. The Fuenmayor system is highly linear for low frequencies up to 0.08 days−1(signals with period greater than 12.5 days). For higher frequencies, the linearity decreases slowly. Compared with others springs analyzed in the literature, Fuenmayor presents a good linearity, similar to Torcal of Antequera (Padilla and Pulido, 1995), and in contrast with Torremolinos (Andreo et al., 2006), Simat (Mangin and Pulido, 1983), Sant Josep (Esteller et al., 1996) and Ancón (Rodríguez et al., 1994). 4.5.1.2 Wavelet Coherence Fig. 4.15 shows the wavelet coherence (Sec. 2.7.2.3) between the rainfall and the discharge at Fuenmayor spring, following the expression 2.44. In each point of the frequency-time map, the coherence result is a value between 0 for dark blue color, and 1 for dark red color. A value 1 for a given frequency band, indicates that the response energy is 100 percent, due to the input rainfall signal enters in a linear relationship. The 5% of significance levels are plotted with black thick line. The left axis is the wavelet scale (hours). The bottom axis is time. Curved lines on either side indicate the COI. The spectral strength spreads out from weak values(blue) to strong ones(dark red). As expected, the wavelet coherence (Sec. 2.7.2.3) shows wide areas of a strong correlation between the input and output signals, with values close to 1, after applying the Eq. 2.44. However, the spectral coherence presents several blue areas (i.e. with very 127 4.5 System identification FUENMAYOR AQUIFER low coherence) in some areas, inside the COI. This suggests there are two behaviors of the system: •Strong coherent behavior that happens when the system clearly reacts in consequence of the rain. •A weak or non-coherent behavior that occurs when it rains after a dry season (0-6,1518, 2830 months). There are periods in the range of June-October 2002, June-October 2003 and MayOctober 2004, in which a non-coherent behavior occurs, showing the non-linearity of the system. Identified empty intervals or dark blue areas correspond to dry seasons. When it rains after a dry season, the soil and epikarst zones need to be wetted, and its maximum water storage replenished before the effective rainfall has some influence on the groundwater and consequently on the discharge output. Fig. 4.15: Coherence spectrum between rainfall and discharge in Fuenmayor spring. The arrows represent the phase angles between the effective rainfall and the discharge in the respective frequency band In wet periods, the linear response of the system and high value of coherence can be easily observed. Conversely, when the system is dry, the reaction to the rain event is not linear. At this point, wavelets can provide their essential support. Thus, the wavelet coherence is useful to determine in which period of time it is recommendable to apply linear methods and when to consider the system as a non-linear model to be analyzed 128 FUENMAYOR AQUIFER 4.5 System identification with different non-linear tools. In conclusion, two main modes have been identified in the system: the wet mode with a linear behavior and the dry mode with non linear behavior. This variation of the coherence evidences the nonlinearities of the system. 4.5.2 Parametric system identification 4.5.2.1 Identification of the impulse response The next step is the identification of the impulse response of the system. Two of the methods presented in Sec. 2.7.1.2 have been applied in Fuenmayor: the Wiener Hopf summation equation (Eq. 2.33) and the Least Squares method (Eq. 2.36). As the estimation is performed by cross correlation, where the estimated impulsional response is near zero after 100 days, only the first 2400 samples (100 days) of the kernel or impulse response have been calculated. It should be stated that both methods presented computational problems with excessive execution times (several hours in an Intel Quad Core at 2.4 GHz) and high memory demands (with 2 GB of physical memory the application filed out of memory). Thus, a decimation process was applied to convert the hourly series into another series with a sample for every twelve hours. After the decimation the calculus was possible. To compare the simulated with the measured discharge, an interpolation process was applied. Fig. 4.16 shows the estimated kernels and the simulated discharge for both methods. The model obtained by the Wiener-Hopf summation equation presents an efficiency (Sec. 2.8.2 ) of E = 0.7941, calculated by Eq. 2.48. The model obtained by error minimization shows E = 0.7934. Thus, both methods are very similar and reproduce the Fuenmayor behavior quite well. 4.5.2.2 Identification of transfer functions Applying to the Fuenmayor spring some parametric identification methods, introduced in Sec. 2.8, and proposed by Ljung (1999) for linear system, such as ARX, OE, ARMAX and BJ, several model structures has been calculated with different orders. The model with the best efficiency (Eq. 2.48) is an OE model with two poles, one zero and seven delays. Applying Eq. 2.24 the following transfer function is obtained : Qs(z) Qu(z)=0.01289 −0.01286z−1 1−1.9882z−1+0.9882z−2z−7(4.1) 129 4.5 System identification FUENMAYOR AQUIFER Fig. 4.16: Kernel estimation of Fuenmayor karst system based on: (a) Wiener-Hopf summation equation and (b) error minimization. As explained in section 2.5.1, poles and zeros are complex quantities that determine the dynamic behavior of the aquifer system. In this case, the TF poles are z1 = 0.9994 and z2 = 0.9888 corresponding to recession coefficients α1= 0.014 days−1and α2= 0.2712 days−1. The first pole value (the slow pole) is similar to the last depletion coefficient, calculated by Trilla and Pascual (1974) using the Maillet method for drought periods. The function also has seven delays (seven units of time). That is the observed lag, around 7 hours, in the response of Fuenmayor to the rainfall. Fig. 4.17 shows the Transfer Function impulse response and the relationship between measured and simulated discharges. According to expression of Nash and Sutcliffe (1970) (Eq. 2.48), the coefficient E evaluate the quality of the simulations obtained for this new model, E = 0.8164. This is the best efficiency of all the models tested. Moreover, it is a manageable model with only five parameters that can be used for analysis, simulation and prediction purposes. A standard procedure in system engineering is to analyze the system model through a set of parameters characteristic of the unit step response. For example: •Static gain = 3.941. Ratio of the system output and the system input values under the steady state condition. •Rise time = 133 days. The time required for the output to reach 90 % of the steady state final value. 130 FUENMAYOR AQUIFER 4.5 System identification Fig. 4.17: Impulse response and simulated runoff in the parametric identification of Fuenmayor spring . •Overshoot = 0. Ratio of the maximum peak value (minus the steady state final value) and the steady state final value. Fuenmayor spring does not show oscillating behavior because both poles are real. The use of these kind of parameter may enable a new comparative analysis of the karst system to be carried out in a standard way. The prediction power of the proposed TF function has been validated with new data set for the year 2010. Fig. 4.18: Prediction of discharge in the first half of year 2010, obtained by parametric identification of the transfer function. 131 4.6 Nonlinear Identification FUENMAYOR AQUIFER A model with an efficiency factor of E = 0.8164 may be valid for many applications in karst hydrology. However, in some cases, a deeper study must be performed to obtain more precise models or analyze the discrepancies with the observed behavior. The identification technique, chosen at this part of the thesis, assumes that the system is linear and time stationary. Fig. 4.18 shows the predicted and measured discharge for that period. 4.6 Nonlinear Identification The prior knowledge on physics of the system, and the previous results from linear identification and spectral analysis, can address the tuning a nonlinear model structure. With one output, i.e. the discharge, and the precipitation as unique input, obtained from Fuenmayor aquifer during 2002-2005 period, three models based on HammersteinWiener have been tried, with different configurations of the component blocks, according with Sec. 2.9.1.2. For the first model M1, the linear block is a linear transfer function corresponding to the orders nb = 4, nf = 2, nk = 1; where nb is the number of zeros plus 1, nf is the number of poles, and nk is the input delay. For the nonlinear blocks, the input and output nonlinearity estimators, two piecewise linear functions with 10 units have been chosen. After performing the test for this model M1, a Nash-Sutcliffe eficiency of E=0.74 was obtained by Eq. 2.48. Asecond model M2 based on the same structure Hammerstein-Wiener was tested. In this model the linear block is an Output-Error polynomial using time domain data, with nb=2 (order of the B polynomial + 1), nf=3 (order of the F polynomial), and nk=1 ( number of samples). The nonlinear blocks, are two piecewise linear functions with 10 units. After implementing the test with M2, the Nash-Sutcliffe eficiency is E=0.81 Third model M3 is also a Hammerstein-Wiener structure. In this model the linear block is a Output-Error polynomial using time domain data, with nb=2 (order of the B polynomial + 1), nf=3 (order of the F polynomial), and nk=1 ( number of samples). The first block is a wavenet and the second one is a piecewise linear functions with 10 units The running test draws out a Nash-Sutcliffe eficiency of E=0.9276 Therefore, model M3 is the best one obtained for the next stage about prediction, in which a good model is required to forecast the behavior of the spring in a long horizon prediction. Fig. 4.19 illustrates graphically the verification of model M3. Moreover, this result leads to those obtained in models of linear identification. 132 COLLINS GLACIER 5.2 Geographical framework 5.2.2 KGI Climate The Antarctic Peninsula is one of the most rapidly warming locations on Earth, with an increasing 0.56◦C per decade measured at the Faraday/Vernadsky station over the last 50 years (Turner et al., 2005). The Antarctic and Subantarctic regions have a prominent climatic interannual variability, compared to other latitudes, due to a number of feedbacks that result from complex interactions between atmospheric circulation, oceans and cryosphere (King and Turner, 1997). KGI climate is determined by the passage of successive cyclonic systems, transporting relatively warm and humid air, and bringing strong winds and heavy precipitation. Ocean is the regulator of the process, so that, in winter, water masses are warmer than the lowest atmospheric layers, thus, the boundary layer warms. In the summer, currents flowing from the Bellingshausen Sea reach the South Shetlands islands and are divided into two branches, one that flows towards the Drake Passage and the smaller one that is addressed to the Bransfield Strait. The low temperatures of these superficial waters decrease the air temperature in KGI (Jiahong et al., 1994). For this reason, although the whole Antarctic continent has an average yearly temperature below zero, the edge of the Peninsula and the surrounding islands are not far from the isotherm at 0 ◦C. Mean annual temperature of permafrost is -1 ◦C at KGI, meanwhile it is between -8 and -10 ◦C at the continent. Southern Annular Mode (SAM) —also named Antarctic Oscillation (AAO)— is an Earth’s atmosphere pattern to describe variability not associated with the seasonal cycle. SAM varies on timescales as fast as weeks. The high and low index polarities reflect the extremes of a normally distributed frequency distribution. SAM is linked to variations in temperatures over Antarctica, sea-surface temperatures throughout the Southern Ocean, and the distribution of sea-ice around the perimeter of Antarctica (Thompson and Wallace, 2000). Therefore, the SAM’s variability affects to KGI’s climate. Ruckamp et al. (2011) investigated the area change and current surface lowering on King George Island, using satellite time series to show changes in glacier extent since 1956 to present for the entire island and differential global positioning system (DGPS) measurements to determine ice surface elevation changes on Bellingshausen Dome. They extrapolated the findings to the entire island in order to provide melt water budgets for individual catchments. Also, Ruckamp et al. (2011) showed a strong interannual variability of two consecutive years by net accumulation measurements and forecasted the disappearance of Bellingshausen Dome by measuring surface lowering rates (1997/98 to 2010) and the ice geometry. South Shetlands the air temperature is characterized by a high value of variability from year to year. A distinct 5.3 years cyclicity can be found in the occurring temperatures. The 139 5.3 Collins glacier COLLINS GLACIER mean annual air temperature meassured by Kejnak (1999) over the period 1944 -1996, with data from the Deception and Bellingshausen stations, was -2.8’C. Its value varied from 0.8’ C in 1989 to -5.2’ C in 1959. The highest temperatures on the King George Island occur in the Admiralty Bay region. 5.3 Collins glacier Glacial geology of Collins comprises the study of the landforms and sediments created by the glaciers, both past and present. The study area covers the south side of the ice cap Collins, commonly called dome Bellinsghausen, located in the southwest of King George Island, and near the Uruguayan Artigas Antarctic Base. Collins glacier, with 1313 km2of extension, occupies almost all of King George Island, except by the south-western end, where the Fildes Peninsula is located. Collins glacier was mapped by Australian National Antarctic Research Expeditions (ANARE) from air photos taken in 1956 and 1960. Collins glacier is named by Antarctic Names Committee of Australia (ANCA) for N.J. Collins, senior diesel mechanic at Mawson Station in 1960. All analysis is carried out on the registered data obtained from at the station of the Insular Antarctica installed by Eraso and Domínguez (2001a). The focus of this work is on Collins glacier data related to discharge and air temperature, although other time series have been also studied, as humidity, pressure, radiation, and precipitation. The data were obtained from 2001 to 2011, at monitoring station of the Insular Antarctica CPE-KG-62◦S by Eraso and Domínguez (2001b). 5.3.1 Hydrological overview The glacier has several domes from which the ice flows to the sea. One of the smallest is Bellingshausen Dome, located in the west end of the island, has several limbs that melt before reaching the sea, and generates nine streams that flow into both sides of the coast, by the slope of the Drake Straits in North side of the island, and through the slope of Bransfield by South. Discharge into the Bransfield is composed by five streams. After following diverse proglacier routes, the whole flow converges at a pond, where an ionospheric observation station is located. The pond output is a unique stream that flows under the bridge leading to the Uruguayan Base Artigas, and finally it reaches the sea at Maxwell Bay. Since this river is the only responsible of the glacier discharge, under the 140 COLLINS GLACIER 5.3 Collins glacier bridge (62◦11’ 03” S, 58◦54’ 41”W), GLACKMA project team installed an experimental station for monitoring the glacier, named as CPE-KG62S. The surface of the selected glacier catchment area, corresponding to the SW sector of the Small Dome or Bellingshausen Dome, was established on the studies carried out with radio-echo sounder. However, its cartography was completed with GPS readings of the perimeter moraine and the surfacing nunataks in this SW sector of Low or Small Dome (Fig. 5.3). The surface of the glacier catchment area (Fig. 5.1 and 5.3) is about 1.31 km2, the development of the peripheral moraine occupies 0.25 km2, and the fluvial catchment extends throughout 1.36 Km2. Then, the total surface of the glacier is 2,92 km2 (Domínguez and Eraso, 2007). The maximum of the ice thickness (over 300 m) occurs in the western and central part of the King George Island ice cap [Pudelko, 2003]. Fig. 5.3: Topographic sketch of ice thickness at Low dome also known as Bellinghausen dome, Collins glacier. 141 5.3 Collins glacier COLLINS GLACIER 5.3.2 Instruments and facilities 5.3.2.1 Discharge volume measurement The discharge volume measurement from the glacier-fed rivers is very difficult to perform accurately due to the rough conditions and the continuous changes of flow regime. An alternative method is the indirect gauge by a pressure-measuring sensor which supplies level values. The measuring station CPE-KG62S has been equipped with register instruments of SEBA Hydrometrie Company, from Germany. The material used at in the first years to produce the time series, consisted of a multi-parameter pressure-measuring sounder (MPS) with sensors for water temperature, conductivity and level. After continuous measuring with hourly records during two complete years, flaws in the external data-logger appeared due to the hard meteorological conditions during the winter months. This originated a non-valid register during the summer of 2003/04, so the sounder had to be replaced by an MDS type sounder that can only measure the river level, but it is an exceptional sounder regarding its resistance in harsh and extreme conditions. MDS-Insider-2 provided by SEBA is an advanced sensor for registration of water level through a high accuracy pressure transducer, with a precision of 1 mm, in the range from 1 cm to 10 m. MSD sensor has been extended to become a multiparameter probe equipped with: •A water temperature sensor, with a range of -20 ◦C to + 70 ◦C and a precision of 0.01 ◦C. •An electrical conductivity sensor, with a range from 0 to 1 mS/cm and an accuracy value of 0.5 µS/cm. Sensors are protected by a cylindrical steel grating attached to probe body, a camera with a watertightness up to 150 meters deep. Other ancillary material complete the station, as a data-logger to save the time series from probes with 64 kbytes memory, lithium and alkaline batteries, M-BUS Interface for laptop PC connection to download data and reprogramming, cable shielded waterproof outer, and SEBA special software for the first inspection of the raw time series generated. 5.3.2.2 Weather measurements About weather measurements, Antarctic environment means a hard challenge for data collection. So much so, King and Turner (1997) suggest that important observations should be complemented by the ability of the human observer. So, the best option is 142 COLLINS GLACIER 5.4 Analysis of time series to take meteorological information from the ground-based station of Bellingshausen (24 m asl). Data can be valid for the Collins study, due to the proximity (about 4 km). The average temperature around the station in August is -6.8 ◦C and +1.1 ◦C in February (Fig. 5.4). From the multiples variables measured in Bellingshausen base, only the following meteorological (input) parameters have been retrieved for the Collins study, though the air temperature will be the prominent parameter to consider in this thesis: •Precipitation (P) in mm. •Atmospheric pressure ( Pa) in mbar. •Global solar radiation ( Rs) in W/m2. •Relative air humidity (Hr) in %. •Air temperature ( Ta) in oC. Air temperature is assumed as the main weather parameter in the Collins glacier study. Fig. 5.4: Bellingshausen monthly temperature statistics: this plot shows the mean, quartiles and range for each month 5.4 Analysis of time series In order to find causality between the weather phenomena and the glacier system, some meteorological variables have initially been considered: precipitation, air temperature, atmospheric pressure and humidity. Each one might present a different contribution to the 143 5.4 Analysis of time series COLLINS GLACIER glacier discharge. After performing XWS (Eq. 2.42) between Collins discharge and each one of the mentioned weather variables, weak correlations regarding discharge were observed between some of them; nevertheless, these extended analysis is left out of this study. The highest correlation was found between the discharge and the air temperature, although closely followed by radiation. The issue about which of the two phenomena have a more direct influence on the melting has been discussed in Sec. 3.3.1. Then, in the next sections it is assumed that the influence of the air temperature on glacier internal ablation is the direct responsible of the glacier discharge. 5.4.1 Time series preparation Data caught from sensors, undergone to extreme conditions of the environment, accidental malfunction of the instruments, or other factors, can supply time series with abnormal values, that can distort the later analysis. For example, some negative values in the discharge time series. The test explained in (Sec. 2.4.2.4) has been applied to time series of this glacier. The result is that outliers have been replaced by regenerated values, according to prediction, trending and smoothing algorithms. 5.4.2 Discharge As mentioned in 5.3.2, at the monitoring station, river levels are gauged to obtain the discharge values. To convert water levels into discharge channel, a customized curve should be designed for this specific glacier, called calibration curve f(H,Q). It is a straightforward solution from the logistic point of view, which associates the river level (H) with the drained volume (Q) (Domínguez and Eraso, 2009). In order to obtain the constants of this function, Domínguez and Eraso (2009) have carried out campaigns of gaugings by fording in the river bed, by covering the greater possible rank levels, with special attention to both minimum and maximum values. The material used was an universal propeller F1 of SEBA company, which is specifically designed to this kind of tasks. Since rivers originated from glaciers presents important variations in flow and load of sediments, before taking each measurement, the river depth should be verified. The suitable place should be a straight stretch with small variations in the cross-sectional profile of its bottom and absence of big obstacles that may distort the speed distribution of the flow. The river width is divided into sections where water speed is measured, which multiplied by section turns out the volume per time unit. The total volume will be the 144 COLLINS GLACIER 5.4 Analysis of time series sum of the volumes of each section. The distribution of the flow speed in a free lamina river decreases exponentially according to the depth, from the surface to the bottom. Domínguez and Eraso (2009) have established the relation between water level and discharge specifically established for the Collins glacier, with the Eq. 5.1. The correlation coefficient R2=0.99 after painstaking works. Qg=0.0159e13.98H(5.1) The discharge measurements carried out in some subpolar stations are often taken with a reading interval of a day, for saving battery consumption, but the loss of information could be too important. For this reason, the values in this Collins multifunction station has been registered in time series with intervals of 1 hour. 5.4.3 Air temperature Because the ice cap is a strong heat sink, vertical and horizontal temperature gradients are constantly varying. For that reason, near-surface temperatures in the glacier are very sensitive to changes in the low layer of the atmosphere. Fig. 5.5 is the representation of the air temperature and their trend. Only in the tenth year, variations in the atmospheric temperatures reaches a half degree. The rate is about one degree per each 20 years. Blue line is the air temperature time series, the red one is the smoothed curve after applying a wavelet filter, and the green one is the linear trending. This is in agreement with King and Harangozo (1998) who estimated that at the last 50 years of twenty century, the warming trend in the west sector of the Antarctic Peninsula, especially between latitudes 65 ◦and 70 ◦S, was approximately 2 ◦C. The same result was obtained by Ferron et al. (2004) in King George Island in the period from 1947 to 1992. Analogously, the trend of discharge has been calculated (Appendix B.5), and the result shows an increasing similar to the air temperature. 5.4.4 The involvement of air temperature in the discharge The time series of discharge and air temperature range from January 21, 2002, until May 31, 2011, with one-hour sampling periods (Fig 5.6). Both time series show a strong yearly cyclic behavior. In cold periods, obviously there are not observable effects in the discharge time series, since there is no flow due to a complete freeze. Therefore, the dynamics is going to be studied only in periods when the glacier is active, i.e. in those 145 5.4 Analysis of time series COLLINS GLACIER Fig. 5.5: Trend of air temperature in the period considered, with a rate of one degree per decade. which the discharge presents a significant reaction to the air temperature. In catching data in general, due to hard conditions of environment or other factors, sensors can provide raw data with some questionable values, i.e. outliers (Sec. 2.2.2). The two series, that are going to be treated in the following sections, require a prior preparing work to check the robustness of the data and perform the proper correctness. The criterion to inquire and revise outliers in air temperature and discharge time series has been the WaveletRosner test (Sec. 2.4.2.4), that also was applied in karst aquifer discharge by Chinarro et al. (2011). The Fig. 5.6 displays the cycle 3 corresponding to the 2013-2004 season hidden by a shaded shape, given that data were lost that year for a sensor breakdown. 5.4.5 Time series power spectrum The graphical output of Fig. 5.7 has been generated by a computer application specifically designed for this analysis, with algorithms based on the work of Grinsted et al. (2004). The basis function for the wavelet transform is a Morlet mother, commonly used in 146 COLLINS GLACIER 5.4 Analysis of time series Fig. 5.6: Measurements of the runoff from glacier Collins at station CPE-KG-62S. a) Time series of the air temperature. b). Idem of the discharge. The third cycle data is lost for failure of sensor. hydrology (Chinarro et al., 2011) because it permits a reasonable precision in situating at timeline, the signal frecuencies (Foufoula-Gergiou and Kumar, 1994; Schaefi et al., 2007). Lafrenière et al. (2003) also uses the Morlet wavelet for analysis of runoff regimes of a glacial watershed. In the ordenate axis, numbers outside the parentheses indicate the equivalent in days. The abscissa axis represents the timeline labeled quarterly, although the calculation of spectral power is held every hour. Values of the spectral power are represented by a distribution of color hue depending on the moment in time series and the frequency (or scale). The cone of influence indicates the useful area of the spectrum without influence of the edge effect. Both spectrograms, discharge (Fig. 5.7) and temperature (Fig. 5.8), show the annual cycles, identified as a horizontal band of high power over 365 days. The discharge also has another high power band corresponding to six months. This band is significantly less powerful in the temperature spectrogram. This is related with the fact that the glacier changes its activity every six months: completely frozen to active and vice-versa. All these results are corroborate with Braithwaite and Olesen (1988). Both spectrograms also present a richer spectrum with frequency components in the band of about 51-60 days. There is another high power band in the range of 21-43 days. The reason for these bands is not well understood, and in a further research based on these facts could throw new findings. 147 5.4 Analysis of time series COLLINS GLACIER Fig. 5.7: Discharge wavelet spectrum at the Collins glacier. The vertical axis has a logarithmic scale in days (hours). The abscissa axis represents the time line labeled quarterly. Values of the spectral power are depicted by a colour map. The cone of influence is also represented to indicate the useful area of the spectrum without influence of the edge effects. 5.4.6 Periods of the annual cycle The structure of discharge time series observed comprises high discharge periods of the glacier that are alternated with winter periods in which low temperature keeps the solid state of the glacier and yields null values in the sensor measurements. The differences observed between the correlograms of the temperature (system input) and the discharge (system output), are also highlighted in the wavelet coherence spectrum. In cold periods, the coherence has no high frequency components because the glacier is completely frozen, although temperature presents diverse frequency components in its power spectrum. Nevertheless, when the glacier is active –in summer– the coherence is very high. Since the cycling behavior has clearly been manifested in the coherence representation, The next stage is to show that a model based on an annual cycle can represent, up to a certain point, any period in the time range considered. Assuming statistically that the probability distribution of a cycle is preserved along different years, a disaggregation study can be posed according with a significant body of hydrological literature (Santos and Salas, 1992). Thus, three periods can be observed in one annual cycle (Fig. 5.9), according to the response of the discharge to the air temperature: 148