Tullio Ricci

Anthony Finizola

Stéphanie Barde-Cabusson

Eric Delcher

Salvatore Alparone

Salvatore Gambino

Vincenzo Milluzzo

Hydrothermal fluid flow disruptions evidenced by subsurface changes in heat transfer modality: The La Fossa cone of Vulcano (Italy) case study

Detecting volcanic unrest is of primary importance for eruption forecasting, especially on volcanoes characterized by highly dangerous, and often seemingly unpredictable, phreatic or phreatomagmatic eruptions. We present a simple and innovative analysis of shallow vertical temperature profiles to depths of 70 cm. These data were recorded at La Fossa cone of Vulcano (Aeolian Islands, Italy), during an episode of increased hydrothermal and seismic activities that occurred between September and December 2009. This work involves the use of the coefficient of determination (R 2 ) on vertical temperature profiles in order to identify changes in conductive versus convective heat transfer modality. The increase in convective heat transfer can be related to the disruption of the hydrothermal system due to its pressurization and/or variation of ground permeability between the hydrothermal system and the surface. While raw temperature data do not evidence any significant variation during the period investigated and the classic temperature gradient is highly influenced by seasonal variations, the fluctuation of R 2 displayed striking spikes that coincided with the seismic swarm inside the volcanic edifice. Such a low-cost device associated with easy real-time data processing could constitute a very promising, yet deceptively simple, technique to monitor hydrothermal systems, in order to assess the hazard posed by high-energy eruptions for populations living close to active volcanoes.

INTRODUCTION

Hydrothermal systems are characterized by important mass and energy transfer through the circulation of hot fluids underground that can be evidenced by geological, geophysical, or geochemical observations. One consequence of hydrothermal circulation is the alteration and weakening of the permeable ground inside volcanic edifices. Hydrothermal alteration greatly decreases the permeability of the medium and disrupts fluid circulation through self-sealing processes. Such modifications may increase pore pressure, the changes in which hydrothermal systems give an indication about the general state of volcanic unrest, and sometimes promote highly hazardous explosive phreatic or phreatomagmatic eruptions [START_REF] Heiken | Phreatic eruption clouds: The activity of La Soufrière de Guadeloupe[END_REF][START_REF] Barberi | A review on phreatic eruptions and their precursors[END_REF][START_REF] Germanovich | The mechanism of phreatic eruptions[END_REF]Ui et al., 2002;[START_REF] Starostin | A transient model for explosive and phreatomagmatic eruptions[END_REF][START_REF] Fournier | Understanding volcano hydrothermal unrest from geodetic observations: Insights from numerical modeling and application to White Island volcano[END_REF].

Studying variations in hydrothermal system activity is crucial to exploring their role in, and potential to forecast highly explosive activity. Therefore, identifying the key parameters that indirectly highlight the state of pressurization of the system is necessary to develop early warning of highly explosive eruptions. Temperature is perhaps the easiest and most obvious observable variable, and can be measured with soil temperature sensors, infrared cameras, or satellite-based multispectral infrared images [START_REF] Gaudin | Heat flux measurement from thermal infrared imagery in low-flux fumarolic zone: Example of the Ty fault (La Soufrière de Guadeloupe)[END_REF].

Unfortunately, the change in temperature in hydrothermal system fluids is not an easy parameter to use in identifying precursors of volcanic unrest. Temperature is affected by many external parameters, such as daily and seasonal variations, rainfall, and wind speed and direction [START_REF] Dawson | Diurnal and seasonal ground temperature variations at Wairakei[END_REF][START_REF] Keszthelyi | Observations of the effect of wind on the cooling of active lava flows[END_REF][START_REF] Chiodini | Carbon dioxide diffuse degassing and estimation of heat release from volcanic and hydrothermal systems[END_REF][START_REF] Hochstein | Measurement of heat flux from steaming ground[END_REF][START_REF] Peltier | Structure of an active volcano associated with a resurgent block inferred from thermal mapping: The Yasur-Yenkahe volcanic complex (Vanuatu)[END_REF].

The aim of this paper is to propose an innovative method to process soil temperature time series that could be used in volcano observatories for real-time monitoring. This new technique can identify disruptions of hydrothermal fluid flow that potentially lead to volcanic unrest.

APPLICATION SITE AND INSTRUMENTATION

The test site where our temperature instrumentation was installed is La Fossa cone on Vulcano Island (Aeolian Islands, Italy), because of its persistent fumarolic activity, the presence of subfumarolic areas (documented in previous studies; e.g., [START_REF] Revil | Inner structure of La Fossa di Vulcano (Vulcano Island, southern Tyrrhenian Sea, Italy) revealed by high-resolution electric resistivity tomography coupled with self-potential, temperature, and CO 2 diffuse degassing measurements[END_REF][START_REF] Barde-Cabusson | New geological insights and structural control on fluid circulation in La Fossa cone (Vulcano, Aeolian Islands, Italy)[END_REF], and the recurrence of fumarolic crises during the past few decades [START_REF] Chiodini | Geochemical and seismological investigations at Vulcano (Aeolian Islands) during 1978-1989[END_REF][START_REF] Granieri | Correlated increase in CO 2 fumarolic content and diffuse emission from la Fossa crater (Vulcano, Italy): Evidence of volcanic unrest or increasing gas release from a stationary deep magma body?[END_REF].

La Fossa cone, formed in the past 6 k.y. on the island of Vulcano, provides the historic basis for the term Vulcanian-type phreatomagmatic explosive eruptions, which are characterized by water-magma interaction [START_REF] Frazzetta | Evolution of the Fossa cone[END_REF][START_REF] Büttner | Identifying magma-water interaction from the surface features of ash particles[END_REF][START_REF] Dellino | Quantitative hazard assessment of phreatomagmatic eruptions at Vulcano (Aeolian Islands, Southern Italy) as obtained by combining stratigraphy, event statistics and physical modeling[END_REF]. The structure of the present-day edifice results from six main phases of activity [START_REF] De Astis | Geological map of the Vulcano Island: Florence, Italy[END_REF], chronologically associated with crater boundaries: (1) Punte Nere, (2) Palizzi, (3) Forgia Vecchia, (4) Pietre Cotte, and (5) Gran Cratere (Fig. 1A). Each of these structural boundaries is associated with relevant temperature anomalies identified through measurements at 30 cm depth [START_REF] Revil | Inner structure of La Fossa di Vulcano (Vulcano Island, southern Tyrrhenian Sea, Italy) revealed by high-resolution electric resistivity tomography coupled with self-potential, temperature, and CO 2 diffuse degassing measurements[END_REF][START_REF] Barde-Cabusson | New geological insights and structural control on fluid circulation in La Fossa cone (Vulcano, Aeolian Islands, Italy)[END_REF], highlighting the major role of these crater ring faults in channeling hydrothermal fluids toward the surface.

Two sites of anomalous temperature were chosen for the installation of our instrumentation that correspond to two structural boundaries located outside the main fumarolic field of La Fossa crater (see Fig. 1A): (1) Gran Cratere (GC in Fig. 1A) and (2) Punte Nere (PN in Fig. 1A), located 250 m and 600 m, respectively, from the last eruptive crater.

We used an Ebro Electronic EBI-2T type 313 data logger to monitor the temperature in the soil. The device temperature range is -40 to +150 °C, with 4 channels and 4 PT1000 sensors (resolution is 0.1 °C and measuring accuracy is 0.2 °C). The temperature sensors were placed at the precise depths of 10, 30, 50, and 70 cm (see Fig. DR1 in the GSA Data Repository1 ) in areas where the soil temperature at 70 cm depth ranged between 70 and 80 °C. It is essential to avoid areas where the temperature signal is saturated by the buffering effect of steam present at boiling temperatures.

Since its last magmatic eruption (A.D. 1888-1890), La Fossa crater has been affected by two main episodes of increased fumarolic activity and temperature. The first occurred during the period 1913-1923 [START_REF] Sicardi | Il recente ciclo dell'attività fumarolica dell'isola di Vulcano[END_REF], and the second began in 1978 after a M 5.5 earthquake [START_REF] Chiodini | Geochemical and seismological investigations at Vulcano (Aeolian Islands) during 1978-1989[END_REF]. Relevant changes in the fumarolic composition occurred in 1979-1981, 1985, 1988, 1996, and in December 2004and 2005[START_REF] Granieri | Correlated increase in CO 2 fumarolic content and diffuse emission from la Fossa crater (Vulcano, Italy): Evidence of volcanic unrest or increasing gas release from a stationary deep magma body?[END_REF][START_REF] Carapezza | Diffuse CO 2 soil degassing and CO 2 and H 2 S air concentration and related hazard at Vulcano Island (Aeolian arc, Italy[END_REF].

At Vulcano, seismicity is associated with volcano-tectonic sources or is related to fluid dynamics within the hydrothermal system [START_REF] Aubert | Double origin of hydrothermal convective flux variations in the Fossa of Vulcano (Italy)[END_REF][START_REF] Cannata | Multiparametric approach in investigating hydrothermal systems: The case of study of Vulcano (Aeolian Islands, Italy)[END_REF][START_REF] Madonia | Thermal anomalies in fumaroles at Vulcano island (Italy) and their relationship with seismic activity[END_REF]. The first case is generally related to the north-northwest-south-southeast Tindari-Letojanni regional fault system dynamics [START_REF] Gambino | Relocation and focal mechanisms of earthquakes in the south-central sector of the Aeolian Archipelago: New structural and volcanological insights[END_REF] and results in a few events per year, while the second comprises seismovolcanic events occurring in the hydrothermal system underlying La Fossa cone, from several hundred to thousands per year [START_REF] Alparone | Time-space variation of volcano-seismic events at La Fossa (Vulcano, Aeolian Islands, Italy): New insights into seismic sources in a hydrothermal system[END_REF][START_REF] Milluzzo | Tornillos at Vulcano: Clues to the dynamics of the hydrothermal system[END_REF]; see the Data Repository).

We compared seismic data with our results, because seismic swarms suggest the disruption of the hydrothermal system, an increase of the permeability inside the edifice, and consequently an increase of hot fluids flow toward the surface [START_REF] Cannata | Multiparametric approach in investigating hydrothermal systems: The case of study of Vulcano (Aeolian Islands, Italy)[END_REF][START_REF] Milluzzo | Tornillos at Vulcano: Clues to the dynamics of the hydrothermal system[END_REF].

During our experiment, seismicity detected by the Istituto Nazionale di Geofisica e Vulcanologia seismic network (Fig. 1B) was characterized by a seismic swarm that occurred at shallow depth below La Fossa cone between 29 September and 16 December 2009 (Fig. 2E). During this seismic swarm, 3471 seismic events were detected in 79 days (an average of 43.9 events/day versus 11.2 events/day characterizing the period before the swarm). Due to the very low energy released, the identification of hypocenters was possible only for the 72 (2.1%) most energetic seismic events, mainly located between 600 and 1200 m below sea level. In terms of number of events per day, this seismic swarm was the most important of the past two decades [START_REF] Harris | Vent temperature trends at the Vulcano Fossa fumarole field: The role of permeability[END_REF]. At the same time, several geochemical parameters increased by one order of magnitude, such as soil CO 2 flux and SO 2 in the plume, indicating disruptions in the underlying magmatic system [START_REF] Inguaggiato | Long-time variation of soil CO 2 fluxes at the summit crater of Vulcano (Italy)[END_REF].

RESULTS

Raw data for soil temperature recorded at Gran Cratere (GC) and Punte Nere (PN) in the period 12 May 2009-28 July 2010 (Figs. 2A and2C) displayed, at both sites, a classic pattern of environmental temperature variations (periodic long-term seasonal variations, and short-term daily variations) disrupted by rainfall events. At the GC site, no temperature variation could be associated with the 2009 seismic swarm. At the PN site, however, only a weak variation of the signal coincided with changes in the seismicity, although there was no clear evidence to support this. In contrast, the temperature gradients calculated between 70 and 10 cm depths displayed a strong seasonality at both sites, preventing any possibility of clearly detecting disruptions related to hydrothermal activity (Figs. 2B and2D). Despite the fact that we can note changes in the temperature gradient coinciding with variations of coefficient of determination, R 2 , at both sites, other changes of the same amplitude in the temperature gradient also appear outside the seismic swarm period. This means that changes in temperature gradient are not exclusively related to hydrothermal system disruptions, and they remain unclear due to a low signal-to-noise ratio.

In order to determine if the 2009 seismic swarm influenced the heat transfer modalities inside La Fossa cone, we computed a linear regression of temperature versus depth and R 2 at each site and time interval. The aim of this analysis was to discriminate the conductive-dominated heat transfer from a more convective one (see the Data Repository); the former is characterized by a linear gradient of R 2 ~ 1, and the latter is characterized by a nonlinear gradient giving lower R 2 values.

Applying an R 2 analysis to our vertical temperature profiles (Figs. 2A and 2C), we obtained the results shown in Figures 2B and2D. At both GC and PN sites, significant anomalous spikes of R 2 coincide with the seismic swarm, showing an obvious disruption in the hydrothermal fluid flow.

DISCUSSION

The comparison between the raw temperature data and the R 2 variations obtained in the same data set displays striking results because no temperature anomaly is obvious in the raw data during the period of seismic swarm. The main external parameters influencing the superficial soil temperature were daily and seasonal temperature variations and contrast between sunshine and cloudy periods. Each one of these parameters is related to pure heat conduction in the soil, and therefore exhibit R 2 ~ 1 despite substantial changes in temperature. These three external parameters have a different influence on temperature values measured by each individual sensor, as shown by the temperature gradient variations (Figs. 2B and2D). Conversely, they do not significantly affect the heat transfer modality (conductive versus convective) of the heat toward the surface, driven mainly by the energy released at depth by the hydrothermal system. Although these external processes strongly modify the temperature values, they do not disrupt the linearity of the temperature gradient as a function of depth. In this sense, the R 2 analysis constitutes an excellent filter for these external parameters and allows for extracting variations related to convective heat transfer of deep origin.

Rainfall is among the external parameters that influence temperature values. Rainfall corresponds to an injection of a cold fluid into the soil, which decreases the temperature and makes the top of the condensation zone deeper (see the Data Repository). As a consequence, the thickness of the conductive zone located above the condensation level increases, as well as R 2 . It can take from several days to a few weeks before a temperature disruption induced by a rainfall vanishes and the temperatures return to their initial values [START_REF] Chiodini | Carbon dioxide diffuse degassing and estimation of heat release from volcanic and hydrothermal systems[END_REF][START_REF] Peltier | Structure of an active volcano associated with a resurgent block inferred from thermal mapping: The Yasur-Yenkahe volcanic complex (Vanuatu)[END_REF]. The duration of the disruption depends on the intensity of the rain and the geothermal energy released from the soil. Although rain events appear to have more significant consequences on the convective component, it is interesting to note that the external parameters previously described (seasonal and diurnal variations and strong sunshine periods) do not disrupt the signal toward a more convective end member. Therefore, a decrease of R 2 is basically the result of an increase in the advective hydrothermal fluid flow toward the surface.

Moreover, unlike remote techniques used for temperature detection (e.g., infrared camera or infrared channels of satellite images), cloudy weather has no adverse effect on the acquisition of soil temperature data, representing a real advantage in thermal monitoring of active volcanoes.

The seismic swarm began 9 days before the decrease of R 2 values, on 29 September and 8 October, respectively. This time delay suggests that the disruption of the hydrothermal fluid flow at the surface could be a consequence of the seismic swarm. Small perturbations temporarily increasing the pore pressure can modify the effective normal stress and trigger seismicity. This phenomenon and other fracturing processes (due to alteration of rock to secondary minerals by reducing the shear stress required to initiate fracturing and/or increases in temperatures in the rock and consequent rock fracturing) are also able to enhance the rock's vertical permeability, thus favoring the rise of fluids (Cannata et al., 2012, and references therein).

CONCLUSIONS

Our results indicate that monitoring the linearity of the vertical temperature profiles at shallow depths (to 70 cm), and therefore the changes from conductive to convective heat transfer in the near surface, is an efficient tool for the identification of volcanic unrest associated with disruption of a hydrothermal system. The installation of such a simple temperature monitoring device could be considered by volcano observatories because it gives fundamental information on changes of heat transfer modality not highlighted by ordinary thermal monitoring. Moreover, data processing can be easily integrated into a real-time monitoring for surveillance purposes. An increase of the hydrothermal activity can be related to pressure and/or permeability changes inside the edifice, potentially leading to explosive activity. Such information is of prime importance for eruption forecasting and hazard assessment and, ultimately, volcanic risk mitigation.

Figure 1 .

 1 Figure 1. Locations of crater boundaries, temperature data loggers, seismic network, and hypocenters of the 2009 seismic crisis at La Fossa cone (Vulcano, Aeolian Islands, Italy). A: Location of the two temperature data loggers over the soil temperature (T ) map (measured at 30 cm depth). Pa-Palizzi, FV-Forgia Vecchia, PC-Pietre Cotte, IVCR-location of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) permanent seismic station at La Fossa cone. B: Location of INGV permanent seismic stations (triangles) operating at Vulcano Island, used to locate the hypocenters of the 2009 seismic crisis, and the seismic events (squares) during the 2009 seismic crisis at La Fossa cone (b.s.l.-below sea level).

Figure 2 .

 2 Figure 2. Temperature variations, coefficient of determination (R 2 ), and seismic events. A: Raw temperature data measured at 10, 30, 50, and 70 cm depth at Gran Cratere (Vulcano, Aeolian Islands, Italy) during the period 12 May 2009-18 July 2010. B: Coefficient of determination (in red) and temperature gradient between 70 and 10 cm depth (in black and yellow) calculated on the temperature data set for the same period for Gran Cratere. z is depth of temperature sensors. C: Temperature data measured at Punte Nere (Vulcano, Aeolian Islands, Italy). D: Coefficient of determination (in red) and temperature gradient between 70 and 10 cm depth (in black and yellow) calculated on the temperature data set for the same period for Punte Nere. E: Number of seismic events recorded at the IVCR station, with the relative cumulative curves. Green area represents the period of seismic crisis evidenced by the cumulative curve of seismicity. The coefficient of determination allows distinguishing the evolution from a conductive to a convective heat transfer induced by a disruption of the hydrothermal system during the seismic crisis.

  l.)

	38.412	100 90	Latitude	38.402 38.412 38.422		A	38.402 38.412 38.422 38.432
				38.392	14.933	14.953	14.973	14.993	38.392
	38.402
	14.953	14.973			14.933	14.953	14.973	14.993
						Longitude

GSA Data Repository item

2015319, temperature monitoring stations, seismicity related to fluid dynamics within the hydrothermal system of La Fossa cone of Vulcano, and heat transfer at shallow depth, is available online at www .geosociety .org /pubs/ft2015.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.

ACKNOWLEDGMENTS

This work has been supported by the University of La Réunion-Institut de Physique du Globe de Paris (IPGP) and the CNRS-INSU (Centre National de la Recherche Scientifique, Institut National des Sciences de l'Univers) project, and a one-month IPGP (Institut de Physique du Globe de Paris) Invited Researcher program for Ricci in the Laboratoire GéoSciences Réunion. Barde-Cabusson was funded by the Spanish National Research Council JAE-Doc (Junta para la Ampliación de Estudios) Program (JAEDoc_09_01319). We thank B. Suski and M.S. Davis for constructive comments concerning the manuscript, and U. Sistometro for field assistance. Comments from reviewers N. Fournier, D. Hasterok, and G. Chiodini greatly improved the manuscript. This is IPGP contribution 3524. Data are available from the authors upon request.


Eli E"'

M Arruda

Mary C Boyer

A THREE-DIMENSIONAL CONSTITUTIVE MODEL FOR THE LARGE STRETCH BEHAVIOR OF RUBBER ELASTIC MATERIALS

A coSSliTt TIVE model i, proposed for tho: deformation of rubber materials which is shnwn tn rcpres..:nt successfully the response of these materials in uniaxial extension. biaxial extension. uniaxial compression. plane strain compression and pure shear. The developed constitutive relation is based on an eight chain representation of the underlying macromokcular network structure of the rubber and the non-Gaussian behavior of tho: individual chains in the proposed network. The eight chain model accurately captun:s the coorx:rativc nature of network d..:formation while requiring only two material parameters. an initial modulus and a limiting chain extensibility. Since these two parameters arc mechanistically linked tn the physics of molecular chain orientation involved in the deformatinn nf rubber, the proposed model rcprcs..:nts a simple and accurate constitUtive model ol rubber ddimnatwn. The chain extension in this network model reduces to a function of the root-mean-square of tho: principal applied stretches as a result of effectively sampling eight orientations of principal stretch space. The results of the proposed eight chain model as well as thos..: of several prominent models arc compared with cxperim..:ntal data of TRELOAR ( 1944. Trans. hm1day Soc. 40, 59) illustrating the superiority. simplicity and pr..:dictivc ability of tho: proposed modd. Additionally. a new set of experiments which captures the state of deformation dependence of rubber is described and conduct~'<.! on three rubber materials. The eight chain model is found to model and prL-dict accurately {he behavior of the three tested materials further conlirming its superiority and ctf..:ctiveness over earlier models.

I. I"'TRODUCTION

THE AIM OF THIS WORK has been to develop a fully three-dimensional. mechanistically motivated constitutive relationship for nonlinear elasticity which would successfully account for the state of deformation dependent response of rubber materials. The need for such a relationship is readily evident to any processor of rubber components or purveyor of rubber goods concerned with aspects of material behavior under large stretch deformation states which arc not simple tension. A good constitutive model should represent the three-dimensional nature of the stress •stretch behavior using a minimal number of parameters to represent physically the deformation process. ldcany, .the parameters should be obtainable from a small number. preferably one, of experiments. Previous models of rubber elasticity such as WA]';G and Gun! ( 1952) and FLORY and REHNRR ( 1943) have described or otherwise accounted for the characteristic 'S'-shaped load versus stretch curve exhibited by rubber materials in uniaxial tension. The load-stretch curve is highly dependent upon the state of imposed deformation, shown, for example, in the data of TRELOAR ( \944) depicted in Fig. I. The authors know of no existing model which accurately represents the behavior of such materials in various deformation states and satisfies the criterion of requiring only a small number of "physically based" parameters. Such material parameters or constants should be independent of deformation state in order to provide a predictive capability to the constitutive model.

The first statistical mechanics approach to describing the force on a deforming polymeric network assumed Gaussian statistics to apply, that is the chains never approached their fully extended length rL = IN where N is the number of statistical links of length I in the chain between chemical crosslinks. [See TRELOAR (1975) for instance for a more detailed description of Gaussian statistics and the corresponding assumptions.] Gaussian statistics yields W= ~G(),f+Jci+lcj-3), G = nk0 (I)

for the strain energy where the lei are the applied stretches and the rubbery modulus G is a function of the chain density, n, Boltzmann's constant, k and temperature 0. Investigators since have built networks from chains described by Gaussian statistics or have modified the chain statistics to allow larger stretches than are afforded by the assumption of Gaussian statistics, then incorporated these non-Gaussian chains into networks of three, four or an infinite number of chains. See WANG and GUTH (1952), FLORY and REHNER (1943) and TRELOAR (1946) for the three chain non-Gaussian, four chain Gaussian and four chain non-Gaussian networks, respectively, or TRELOAR (1975) for an overview of network models including his own proposed method of averaging the contribution of a single chain over a large number of orientations.t These models have in common two physically based parameters, the rubbery modulus G given in (I) and a chain locking stretch A.L defined as the value of the chain stretch when the chain length reaches its fully extended state. Chain stretch is always given by the current chain length divided by the initial chain length fchain

)•chain = --• ro (2)
The initial chain length is given from random walk statistics as r 0 =for Thus for the fully extended or locking chain length, rL = IN, the locking stretch becomes

(3)

Other statistical models have investigated affine versus phantom deformation of junctions in networks and these have been discussed by MARK and ERMAN (1988). In many cases models such as these require a third or fourth parameter to describe some measure of chain interaction in the network.

Excellent historical perspectives of the phenomenological invariant-based models of rubber deformation such as those of MooNEY ( 1940), RIYLIN (1948), V ALANIS and LANDEL (1967) and OGDEN (1972) have been given by TRELOAR (1975,1976). The theories of this genre are aimed at obtaining an expression for the elastic strain energy which is of the form dictated by continuum mechanics as concerns an initially isotropic, incompressible, hyperelastic solid. The phenomenological theories lack a direct physical connection to the underlying mechanisms of deformation. The most sophisticated of these is that of OGDEN ( 1972) which proposes the following form for the strain energy function with fln and the 'Y." as experimentally fitted constants. Ogden makes no attempt to link the adjustable parameters to any physical deformation mechanism. The Ogden model is essentially empirical and generally requires more than one experiment to obtain the number of constants required to capture state of deformation dependence.

tThe averaging was later published in TRELOAR and RID!t'G (1979).

PROPOSED MODEL

Requirements of a model

The existing models referenced above have merit and are frequently used in treating rubber deformation problems. Many of the early models capture the 'S' -shaped load versus stretch curve exhibited by rubber materials in uniaxial tension. A relative dearth of experimental evidence exists in the literature for the response of rubbers in other deformation states up to very large stretches save for the work of TRELOAR ( 1944), in uniaxial extension, biaxial extension and shear, considered by many to be the quintessential rubber data,t and the extensive biaxial experiments of JoNES and TRELOAR (1975). Consequently, many of the existing models fail in the task of describing the response of a rubber material under different states of deformation without changing the model parameters. In other cases the models are prohibitive in the number of parameters required to fit the data or the mathematical complexity they represent. The need for a constitutive relationship which possesses mathematical simplicity, requires one test to characterize the material and has a limited number of parameters, has prompted the authors' development of the following model.

Langevin chain statistics

The statistical mechanics approach to rubber elasticity models the rubber chain segment between chemical crosslinks as a number N of rigid links of equal length /. The rigid link is that segment length of the actual chain which undergoes rigid body motion in response to an imposed strain, depending on the actual rubber material the statistical rigid link may span one or several repeat molecular units. The initial chain length is taken from a random walk consideration of N steps oflength /, and is denoted by r 0 , ro = J!vz.

(5)

The fully extended chain has approximate length IN so that the limiting extensibility (or chain locking stretch), defined as the final length divided by initial length, is given in terms of the statistical parameters as )~L = "/N. At any value of chain length the most probable angular distribution of rigid links about the chain vector length may be found. Following the use of Langevin statistics by KuHN and GRON (1942), the probability of the most probable link angle distribution about a given vector length is taken to be equal to the probability of the vector length. The chain vector length is denoted by rchain• In this way Kuhn and Gri.in obtain an expression for the probability density function for chain lengths and subsequently the configurational entropy of a stretched chain of current length rchain i" Some debate has persisted over whether the characteristic upturn in the observed stress-stretch behavior is solely due to crystallization and therefore. not present in all rubbers. In the following section we offer additional experimental results on three rubber materials which exhibit the same basic. characteristic behavior as that of Treloar's data. ,

W = nk0N --/]+In--. ---ec, Nl smh f3 (7)
where n is the chain density, c' is a combination of constants and the remaining terms are as defined previously. The material is considered to be incompressible and the principal stresses can be determined from the work of deformation to within an arbitrary pressure, p,

(8)
where the A., are the principal stretches and the pressure may be determined from the boundary conditions. The stress-stretch relations are frequently written in terms of the difference in two principal stresses in order to eliminate the pressure term.

(9)

Use of Langevin chain networks

The use of Langevin chain statistics in networks has been considered by [START_REF] Wang | REFERE~CES Anisotropy and localization of plastic deformation[END_REF] in their three chain model by TRELOAR (1946) in his extension of the tetrahedron model of FLORY and REHNER (1943). TRELOAR (1975) andlater TRELOAR andRIDING (1979) also considered a model of a large assembly of chains which were averaged by integrating over many spatial orientations. Analytical results of non-Gaussian three chain, tetrahedron and total assembly of chains models are presented in TRELOAR (1975) for uniaxial tension. The three chain and tetrahedron models have been considered by the present authors for use as constitutive models of rubber deformation, results of simulations of different deformation states using both models will be presented later and compared to the proposed model. The three chain and tetrahedron model systems are sketched in Figs 2 and 3 respectively for undeformed, uniaxial extension and biaxial extension load geometries.

These models were also previously considered for their general predictive capabilities of the characteristic 'S' -shaped uniaxial response by TRELOAR (1975). He observed that the response of the tetrahedron model was slightly dependent upon the orientation of the volume element with respect to the extension direction. In their 1952 paper, Wang and Guth stipulated that the faces of the three chain cube element be aligned with the principal stretch space during deformation. TRELOAR (1975) concluded that these models possess an inherent anisotropy which could be eliminated by a system of a large assembly of chains, each of length equal to the root-mean-square chain length, r 0 = .,jNz, which occupy random orientations in space. In this manner the contributions from chains in several orientations may be averaged to give the overall response. His proposed total assembly of chains model requires integration over the distribution of chains for each subjected stretch state, a procedure which was mathematically prohibitive 15 years ago and remains cumbersome today.

The authors agree with Treloar's assessments of the existing network models except for his reluctance to restrict the orientation of a model to some specific relationship with respect to principal stretch space as Wang and Guth had done. Motivated by Treloar's assertion to present a model which averages several spatial orientations, but recognizant of the need for a mathematically concise representation, we have proposed a model which possesses the cubic symmetry of principal stretch space as it averages eight orientations of that space in determining the network response. The clear advantage of this technique is its ability to simulate a true network response of cooperative chain stretching which presents a clear picture of the deformation process with mathematical ease.

Proposed model geometry

The proposed model considers eight orientations of chains in space which may be envisioned by the eight chain network system sketched in Fig. 4 uniaxial extension and biaxial extension loaded configurations. The chain length, rchain• is sought for each of the chains in the eight chain model as a function of the imposed principal global stretches A~> A 2 , A 3 in order that ( 7) may be differentiated to obtain the stress-stretch relations (9). Consider the cube containing the eight chain configuration in Fig. 5. The cube edges are taken to remain aligned with principal stretch space during deformation, chains linked at the center of the cube extend to the eight corners. In principal space the cube is allowed expansion along each principal direction subject only to incompressibility which may be expressed as (10) Note that the restriction to principal stretch space does not limit the usefulness of the model to axisymmetric deformations. For example, in pure shear the principal stretches are }. 1 = ),, A 2 = 1 and A 3 = A-1 • The directions in which these stretches act rotate continuously with the deformation and this rotation is monitored by the standard kinematics of finite strain deformation making use of the polar decomposition and subsequent extraction of the stretch tensor eigenvalues from the general (1976) or FARDSHISHEH and0NAT (1972). The eight chain model is always stretched in the principal frame and is in general oriented differently from the laboratory frame. The rationale for use of the model in this manner follows from the consideration that in response to any deformation a principal stretch frame exists and the chains in that reference frame will undergo stretches describable by the principal values of stretch, i. 1 • ) 2 and i. )3

(II)
In Fig. 6 the cube is stretched by) J, Jc 2 ] and i.i( so that the cube edges measure }. 1 a 0 ,

). 2 a 0 and ). 3 a 0 in the[, j and k directions, respectively. A chain vector from the center of the cube to a corner may be written for one chain as

(\ 2)
This chain has vector length

( \3)
as do all remaining chains in the given network geometry, regardless of deformation state. Substitution from ( 11) and ( 5) into (13) gives the chain length in terms of the statistical parameters and the principal stretches In this form the expression for the chain vector length is suitable for substitution into (7) which with (9) yields the following stress-stretch relation for the proposed model

nk0 -I [Achain] (Jc I-Jci) (JI -(J? = --N2? -- - 3 jN Achain ' ( 15
)
where the chain stretch, Achain = rchain/r 0 , is the same in each chain and is given by

) 1 ('2 12 )2)1/2 •chain=~ AI +A2 + •3 • • J3 (16)
Because each chain in the system undergoes a stretch equivalent to that in every other network chain, the model is likened to averaging the contributions of a single chain over eight spatial orientations.

The three stretch invariants are given by

II= ).f+A~+).~, I2 = JciJc~+Jc~Jc~+2iJcL I3 = JctJciJc~. ( 17
) (18) (19)
The expression for chain stretch is seen to reduce to a function of the first stretch invariant, I~> and ( 16) may be rewritten as

-I 1/2 Achain -J3 I 1 • (20)
The strain energy of the proposed model may be found from integration of (15) using the series expansion form for the inverse Langevin function given for example in TRELOAR (t 954). The first five terms for the strain energy of this model are (21)

Note that the strain energy in (21) exhibits a nonlinear I 1 dependence as a result of the chain stretch being defined in terms of / 1 only.

EXPERIMENTAL PROCEDURE

The uniaxial compression test and plane strain compression test represent near extremes in the behavior of polymeric networks under biaxial deformation states. FIG. 7. Sketch of plane strain compression testing apparatus. Specimen (A) is compressed by punch (B) in the Y-direction and expands in the X-direction. The die (C) constrains the specimen along the Z-axis.

During uniaxial compression the polymer chains in the material stretch freely in all directions within a plane perpendicular to the load axis. Under plane strain compression the specimen is held within a channel die which completely constrains the material in one direction while allowing expansion in a second direction perpendicular to the first. The material flows along the expansion direction as a result of the compressive load applied perpendicular to the plane containing the constrained and expansion directions. In the plane strain compression die apparatus sketched in Fig. 7, expansion occurs along the X-axis, the material is constrained along the Z-axis and Y is the load direction. Polymer network stretching occurs only along the axis of free expansion in plane strain. The resulting stress versus stretch curve in plane strain compression differs from that in uniaxial compression because of the vastly different chain stretch behaviors in these two compression geometries.

Three commonly available rubber materials were chosen for this study: silicone rubber (40 durometer), neoprene rubber (60 durometer) and gum rubber (hardness unknown).t Uniaxial cubes of length 12 mm and plane strain specimens measuring 12 mm x 9 mm x 12 mm were cut from these materials. rn the plane strain compression tests the 9 mm dimension was constrained. Uniform deformations of the order of -200% strain were achieved in uniaxial compression and -150% strain in plane strain compression. The specimen (and die, in the case of plane strain) was placed between flat plates on the ends of the load strain on an T nstron model 1350 servohydraulic testing machine. Liberal lubrication was applied on all specimen surfaces in contact with plates or die walls to insure uniform deformations; a MoS 2 based paste made by Dow Corning was deemed sufficient for lubricating these materials. Linear ramp tests were performed on all three materials for both deformation states at displacement rates of approximately I mm s •• 1

• The load and displacement output responses of the testing machine were converted to digital signals by a Kiethly AID t These three materials were obtained at the Greene Rubber Company of Cambridge. MA. 4. RESULTS

Explanation of the experimental data

As shown earlier in Fig. 1, the data on Treloar's vulcanized rubber include tests in uniaxial extension, biaxial extension and pure shear to very large extension ratios and are plotted as load normalized by initial area versus stretch. The strong state of deformation dependence of the large and limiting stretch behavior of the polymeric network materials is shown in the data of Figs 1, 8, 9 and 10 and remains a primary obstacle to true characterization of rubber elasticity. Additionally, ARRUDA and BoYCE (1991) have also reported a similar state of deformation dependent strain hardening response in glassy polymeric materials. These data on several different polymeric materials, both above and below the glass transition temperature, clearly demonstrate a certain level of generality to the presence of limiting stretch states in polymeric materials whether rubbery or glassy.

The equi-biaxial extension data in Fig. I varies significantly from either the uniaxial extension or pure shear data. The divergence of the biaxial data is due to the nature of the molecular chain stretching in equi-biaxia\ extension versus uniaxial extension. Biaxial extension offers two venues of principal tensile stretch (}. 1 =A, l 2 = i_,

)., = ). • 2 ); an initially isotropic network of chains will reach limiting chain extension due to stretching in both directions providing a planar state of orientation. In uniaxial tensile deformation the chains extend along one direction only (A 1 = )., }. 2 = ). 12 , ). 3 = ).-12 ), additional stretch is thus allotted through the drawing of material from the transverse directions, and the onset of limiting chain extension is delayed with respect to the biaxial deformation locking. The pure shear data are plotted in terms of the maximum principal stretch, A 1 , vs the corresponding normal force (i.e. not the shear force) which acts in the direction of /. 1 • Pure shear deformation is more closely related to uniaxial extension than to equi-biaxial extension because chain stretch occurs due to stretching along one principal direction (A 1 = },, A 2 = 1, A 3 = }_-1 ) with chains being drawn from one direction transverse to the extension direction in pure shear. Thus a pure shear experiment yields a limiting stretch value which is similar to that obtained in uniaxial extension. A good model should be able to capture the observed differences between biaxial and uniaxial tensile behavior as well as the similarities in behavior between uniaxial and shear tests.

The plane strain compression and uniaxial compression data in Figs 8, 9 and 10 diverge for reasons similar to those used to explain the biaxial extension divergence from the uniaxial extension and pure shear data above. Uniaxial compression offers two directions of equal principal stretch (..1.

1 = A, A 2 = ), • 1 ' 2 , A 3 = )_-1 '
2 ), whereas in plane strain compression chains extend along one direction only. The plane strain compressionstretchstateisequivalenttothatofpureshear(A-1 = ),,A 2 = l,A-3 =A -1 ).

Explanation of simulations

Equation ( 15) has been used together with the appropriate boundary conditions to simulate the uniaxial extension, biaxial extension and shear reponses of the proposed model. Characterization of a particular material requires determination of two parameters, the initial modulus and limiting extensibility, which is accomplished by fitting the model to the data for any one deformation state. The authors chose to fit Treloar's uniaxial extension data which yielded the following parameters for this material: CR = 0.09 and N = 26.5. These parameters were used in the simulations of biaxial extension and pure shear tests. The results appear in Fig. 11. Using constants determined from this single set of uniaxial extension data the model quite accurately predicts the biaxial extension and pure shear data. These results demonstrate the predictive and fully three-dimensional aspects of our proposed eight chain model.

This exercise was repeated for both the three chain model of WANG and GuTH (19 52) and the non-Gaussian tetrahedron model as modified by TRELOAR (1946TRELOAR ( , 1954)). The stress-stretch relations for the Wang and Guth model are of the type

nkEJ j -(. .. 1 { A 1 } • _ 1 { ),2 }) a1 -(52= 3-N A12" J!V -A 2 5£ fo , ( 22
)
where the parameters have the same meaning as in ( 15). The numerical treatment used for the tetrahedron model followed the early method of TRELOAR ( 1954) which allowed the central junction point to seek an equilibrium position for affine displacements of the tetrahedron corners. Equilibrium was found using an iterative Newton scheme after each stretch increment for the condition of no net force on the junction point. Stretch-stress relations are a function of the junction point equilibrium position and must be determined numerically. Each of these models requires one set of data to determine an initial modulus and locking stretch. The constants were found by fitting the models to the uniaxial extension data then used in predicting the biaxial extension and pure shear responses. The results for the three chain model are given in Fig. 12 and for the tetrahedron model in Fig. 13; each figure includes the values used for the model parameters. As these figures clearly show, the early models do not predict much of a state of deformation dependence beyond the initial modulus and cannot be used as constitutive indicators of rubber elastic deformation. The three chain model shows no state of deformation dependence beyond a stretch of 4.0 for the three stretch states examined; the results of the simulations for three states of deformation converge at large stretches whereas the data diverge. The failure of the three chain model to capture the state of deformation dependence lies in the absence of a cooperative nature of network chain deformation in this model except through the imposition of incompressibility. In uniaxial extension, depicted in Fig. 2, the chain parallel to the extension direction stretches until it reaches the locking stretch. The remaining two chain stretches are determined by incompressibility and do not contribute significantly to the force. Under equi-biaxial extension, two coplanar chains extend independently until simultaneously reaching the locking stretch and the third chain stretch is determined by consideration of incompressibility. Therefore the limiting stretch in equi-biaxial extension as predicted by the three chain model has the same value as that obtained in uniaxial extension and in pure shear, the result is convergence of the load vs stretch predictions for the three different stretch states rather than the actual divergence of these curves which the data show.

The tetrahedron model shows a slight state of deformation dependence at small stretches which is retained at large stretches. This model has some success in predicting a state of deformation dependence because the chains in this network do respond cooperatively to the deformation, see Fig. 3. Four chains stretch along the extension direction in uniaxial extension resulting in an overall network limiting tensile stretch which exceeds the locking stretch of a single chain (A.L = J'N). In biaxial extension perpendicular to one chain in the tetrahedron the remaining three chains undergo stretching. The result is again increased network extensibility over the single chain limiting stretch to a value which differs from the tensile locking stretch. The different values for locking in uniaxial and biaxial extension are manifested in the load-stretch predictions of the four chain model which reach asymptotes at different values of stretch for different states of deformation. The tetrahedron model does not possess symmetry with respect to the principal stretch space. The load versus stretch predictions depend upon the orientation of the four chain tetrahedron with respect to i. 1 , i. 2 and ). 3 in addition to the magnitudes of the ..:t,.

The eight chain model is clearly superior to the other statistical models in the state of deformation dependence predicted. The advantages of the eight chain model arc its proximity in kind to a system which is initially isotropic, the use of Langevin chains which capture the effects of limiting chain extensibility and a network configuration which responds cooperatively to an imposed deformation. These features allow for mechanically simulating the state of deformation dependence of the behavior during large deformation while retaining mathematical simplicity. The features of the eight chain model as seen in Fig. 4 combine the favorable aspects of the previous three chain and four chain network models. The model possesses symmetry with respect to the principal stretch space which was required by WANG and GUTH (1952) in their discussion. All eight chains stretch uniaxially in response to an imposed uniaxial extension deformation as is true for the four chains in the FLORY-REHNER (1943) model. Again the limiting network stretch of the eight chain network exceeds the extensibility of a single chain in extension. In biaxial extension all eight chains extend identically due to stretching in two principal directions. The resulting limiting stretch is something other than either the individual chain locking stretch or the limiting stretch in uniaxial extension. A significant feature of the eight chain model is that all chains stretch equally under biaxial extension as well as uniaxial extension. The chain stretch is always the root-mean-square of the global principal stretch state. This unique feature of the eight chain model is true for all other deformation states; all chains stretch equally in response to any given deformation state.

Additional existing models were similarly tested against Treloar's data. These results are discussed in the Appendix. state of deformation dependence because of inabilities of the models to effectively represent a network response. In the case of the eight chain model the successful prediction is a result of correct usage of the finite extensibility parameter as that of an averaged measure of applied stretches over eight spatial orientations. The mathematically tractable eight chain model of ( 15) represents a true constitutive relationship for rubber deformation which is able to correctly account for the large stretch deformation response in one state of deformation, then with that characterization, predict other states of deformation.

A polymer below its glass transition temperature exhibits the same network response in strain hardening as rubber materials do during large stretch deformations. ARRUDA and BoYCE (1991) have shown the eight chain model to be successful in predicting the state of deformation dependence of the strain hardening response of glassy polymers whereas the three chain and tetrahedron models again fail.

SuMMARY

A physically based constitutive model for large stretch rubber deformation has been proposed which has been specifically designed to account for the three-dimensional state of deformation dependence in networked solids. The eight chain model presented here is formulated such that the nature of the state of deformation dependence is clearly seen to be the result of a network of chains reaching the individual chain extensibility limit at different imposed global stretch levels for different stretch states. The eight chain model successfully accounts for the state of deformation dependence using a rubbery modulus and a locking stretch as its only two parameters, both of which can be determined from a single experiment. Present in all statistical models, these parameters are physically linked to the polymeric network and therefore provide a basis for including other aspects of rubber elastic behavior such as temperature dependence, swelling and Mullin's effect. Indeed,-the rubbery modulus term explicitly includes temperature.

Experiments were conducted on three rubber materials to illustrate the state of deformation dependence of these materials and to determine the predictive capability and effectiveness of the proposed model. In addition the model was compared to existing rubber elasticity models in the ability to capture the response of data in the literature (TRELOAR, 1944) in three states of deformation. The eight chain model was shown to be superior in its overall ability to successfully account for the threedimensional nature to the underlying mechanics of network solid deformation for a total of four materials representing a range in material hardness and extensibilities.

The eight chain model contains many attractive features of concern in modelling the complicated deformation procedures involved in the finite straining of amorphous polymeric solids. It retains mathematical feasibility as a two parameter model and also contains the ability to respond to an imposed deformation state in a manner which simulates the actual mechanism governing the state of deformation response of rubber materials. An additional statistical mechanics model of FLORY and ERMAN ( 1982) which accounts for chain interactions has been considered. The elastic strain energy of the network is found from the sum of phantom and constraint contributions w = wrh+ w"

where Wrh of phantom Gaussian chains is 6 .+r:O-gd-e-n:-H""1-.-.61-B-:03:----:A,-1--1-.3-, and where K is a measure of the strengths of the constraints which depends on the relative sizes of free (phantom) fluctuations and actual constrained fluctuations. This model depends on the parameters n, ¢ and K. In it K-> CXJ for completely constrained junctions and K-> 0 in the phantom chain limit, n can take on any large, positive value and ¢ must be greater than two.

Results of simulations with the Flory and Erman model appear in Fig. A2. The model results shown in Figs AI and A2 reveal that of these existing models, only the Ogden model captures the state of deformation dependence of deformation. However, in order to capture this behavior, the Ogden model required six parameters.

  FIG. I. Data from TRELOAR(1944), plotted as force per unit unstrained area versus stretch, showing the state of deformation dependence of a rubber material in uniaxial extension, biaxial extension and shear.

  c-N -iff f3 + ln sinh 7J ' (6) where cis a constant, k is Boltzmann's constant and f3 is the inverse Langevin function, f3 = 2' -1 [rchain/ Nl], for the Langevin function defined as 2'[/3] = coth f3-(1 / /3). The use of Langevin statistics properly accounts for the limiting chain extensibility. The work of deformation is proportional to the entropy change on stretching the chains from the unstretched state and may be written in terms of the chain length as ( r chain f3 )

  FIG. 2. Three chain rubber elasticity model for (a) undeformed, (b) uniaxial extension and (c) biaxial extension configurations.

  FIG. 3. Four chain rubber elasticity model for (a) undeformed, (b) uniaxial extension and (c) biaxial extension configurations.

  FIG. 5. The unstretched network for the proposed eight chain model.

  3 .The unstretched network includes eight chains of length r 0 = }Nl inside a cube of dimension a 0 . From this geometry 2 Clo = --= ro.

  FI<i. 6. The eight chain network in a stretched configuration. rchain = ~jNl(Jcr+Jc~+Jcj) 112 • (14)

FIG. 9 .

 9 FIG. 8. Ramp tests in uniaxial compression and plane strain compression for silicone rubber.

  . Ramp tests in uniaxial compression and plane strain compression for neoprene rubber.

  FIG. 11. Results of simulations using the eight chain model versus data from TRELOAR (1944) in uniaxial extension, biaxial extension and shear.

FIG. 12 .

 12 FIG. 12. Results of simulations using the Wang and Guth model versus data fromTRELOAR (1944) in uniaxial extension, biaxial extension and shear.

  Results of simulations using the four chain model versus data from TRELOAR ( 1944) in uniaxial extension. biaxial extension and shear.

  FIG. 14. Results of simulations using the eight chain model and data on silicone rubber.

  17

  FIG. I 5. Results of simulations using the eight chain model and data on gum rubber.

  FIG. 16. Results of simulations using the eight chain model and data on neoprene rubber.

  FIG. 17. Ditect comparison of the uniaxial compression responses of silicone, gum and neoprene rubbers.

  enological models were considered, the well knownMooney-Rivlin (1948) relation and theOgden (1972) model. The strain energy in the Mooney-Rivlin relation is given as W= C,(l,-3)+C 2 (/ 2 -3), (AI)where C, and C 2 are constants. The strain energy expression for the Ogden model was given in the text of this manuscript. The results of simulations with both of these models appear in Fig. A I along with the constants needed for each model. The simulations with the Ogden model were taken directly from OGDEN (1972). A minimum of six independently adjustable parameters are required by the Ogden model to fit these three deformation states shown in Fig. A I ; more than one deformation state may have been considered in determining the six constants as only the first four constants are necessary to produce the uniaxial extension result.

H2

  FIG. A\. Results of simulations using the Mooney••Rivlin and Ogden models versus data from TRELOAR (1944) in uniaxial extension, biaxial extension and shear.

AcK ~OWLEDGEMENTS This work is sponsored by The National Science Foundation (MSM-8818233) and the MIT Bradley Foundation. ARRUDA, E. M. and BoYCE, M. C. 1991 F ARDSIIISHEH, F. and ONA T, E. T. 19T2 FLORY. P. J. and ERMAN, B.


Hubert Godfroy

Jean-Yves Marion

Abstract Self Modifying Machines

Introduction

Self-modification shifts the current programming language theory paradigm. In textbooks, see for example Winskel's [11], programs are defined by fixed text. The (text of the) program is available and invariant with regard to any execution. The theory leans on a clear distinction between programs and data. This separation between data and programs is also reflected in models of computation. In the case of Turing machines, or of Random Access memory, the program is "wired" inside the control unit. The case of abstract machines for functional languages is a bit different and may be seen as an exception. Indeed, functions are seen as values and data as terms (e.g. Church numerals) but this apparent confusion between data and programs is an artefact rather than a fundamental property of abstract machines. Indeed, there is (always) a clear distinction between data and terms at higher semantics levels.

There are models of computation in which data are first order citizens like Random access stored programs (RASP) [4,5]. However, RASP are not satisfactory because self-modifying behaviors are hidden.

Our objective is to devise a model of computation, which has self-modifying mechanisms. There are at least two dual motivations. The first is to compile a program into a self-modified code. Nowadays, a lot of programs protects themselves against analysis by hiding their behaviors thanks to packers. A packer compresses (and encrypts) a program preserving its semantics. The result is wrapped with a short decompression routine. At runtime, this routine unpacks the code and then it transfers the control to the initial code. Unpacker builds a quite simple, but efficient, self-modifying procedure. We could think of much more involved self-modifying mechanisms. For this, we should make a clear distinction between two notions: programs and codes. Programs will be considered invariant unlike codes which are executable data. Given a programming language L , a program is compiled into a code, which is designed to run on a given model of computation. This model of computation may have features which are not in L semantics like the ability of modifying a code at runtime. Our abstract machine ASM 2 would provide an intermediate stage for compilation where selfmodifying mechanisms could be used.

The other motivation is the analysis of self-modifying codes, because nowadays all malware use this kind of obfuscation [2, 1]. The goal is then to recover from codes the original programs in order to understand its meaning.

This work is also motivated by the work of Jones et al.

[6] on biological systems. In this model, computations are made by a net of blobs where data are programs, with a strong notion of locality. There are models of self-modifying codes like Della Preda et al. [9] where a fixed-point semantics of self-modifying programs is defined, or like [3,7] where a framework developed in the COQ proof assistant is constructed. Finally, self-modifying codes are related to Kleene's second fixed point theorem as it is explained in [8].

Abstract Self Modifying Machines

The salient feature of Abstract Self Modifying Machines (ASM 2 ) is the memory architecture. The memory consists of two zones of registers: (i) the executable zone containing code registers and (ii) the data zone containing data registers. The access to both zones are restricted and are different. We can run every register in the executable zone, but we cannot write on it. We can write on every register in the data zone, but we can not run it. This apparent separation between data and codes is broken by transfer instructions. We have an instruction which transfers a data register into the executable zone, together with an instruction which runs any code register. As a result, we can construct a code inside the data zone. When it is ready, we can transfer it to the executable zone and run it. This is quite similar to a decompress/decryption loop currently used by unpackers. As a result, a register contains a word, which may be decoded as a sequence of instructions when it is put in the executable zone. This is a typical behavior of a self-modifying code.

Formally, a memory is a pair of zones (X, D) where X and D are finite sequences of registers referring to the executable zone and to the data zone. Each register is identified by a natural number. Each register contains a word on some alphabet Σ. Concatenation over Σ * is noted •. Moreover, we have a set of instructions I and an encoding partial function enc :

Σ * → I A state (X, D) | RP | IP consists of a memory (X, D)
, a register pointer RP and an instruction pointer IP . The machine counter is broken into a register pointer RP and an instruction pointer IP . Hence, the current instruction starts from address pointed by IP of the register pointed by RP . At each step, the current instruction is run. The next instruction depends on the instruction kinds. The instruction pointer IP then moves accordingly. So, we write IP = IP + 1 to increment the register IP of the length of an instruction opcode, i.e. to move forward of one instruction, and IP = n to point the n th instruction.

The execution flow may move from one register to another by changing the values of the register pointer RP .

If the register n is in executable zone (resp. data zone), we use X

Transition rules

We will now describe the transition relation. The rules are divided in two categories and are given in the shape : (X, D)

| RP | IP inst -→ (X ′ , D ′ ) | RP ′ | IP ′
where inst is the instruction beginning at the address pointed by IP is the code register X[RP ].

Calculation rules Calculation rules are devised to compute in an environment where data and codes are separated. No executable registers are modified.

-We concatenate the letter a at the end of the data register D[n]:

(Move) (X, D) | RP | IP move a,n -→ (X, D[n ← D[n] • a]) | RP | IP + 1
-We erase the letter a on the top of the data register D[n]:

(Pop) (X, D) | RP | IP pop n -→ (X, D[n ← u]) | RP | IP + 1 where D[n] = a • u -We jump to the instruction m th instruction in the register X[RP ]. For this,
IP is set to m:

(Jump) (X, D) | RP | IP jump m -→ (X, D) | RP | m
-If the value on the top of the register D[n] is a, then we run the next instruction, otherwise we jump to the second instruction:

(Case) (X, D) | RP | IP case an -→ (X, D) | RP | IP + 1 if D[n] = a • u (X, D) | RP | IP case an -→ (X, D) | RP | IP + 2 otherwise
Code motions -We run the code register n starting at the first instruction (i.e. = 0):

(Exec) (X, D) | RP | IP exec n -→ (X, D) | n | 0
-The data register n is moved to the executable zone. As a result, it can be executed thanks to the above rule:

(Activate) (X, D) | RP | IP refl n -→ (X ⊕ D[n], D \ D[n]) | RP | IP + 1
-The code register n is put into the data zone. It cannot be run but it can be modified:

(Inactivate) (X, D) | RP | IP reif n -→ (X \ X[n], D ⊕ X[n]) | RP | IP + 1
Typically, program halts if IP "goes below" the current code. There are other reasons to halt. The register IP may point to a word, which is not the opcode of an instruction. Another way to halt is when a jump is performed with an address outside of the register addresses.

Example

We present here a very simple self-modifying program which decrypts the register C with key in register K (simply a XOR), puts result in register M and then executes it.

We suppose that our alphabet Σ is the set {0, 1, #}. We need four registers I, C, K and M whose last tree ones have already been described, the first is the initialising register, which is executed first. In the beginning, C, K and M belong to D and I to X. The register and instruction pointers are initialized to I and α 0 . The first state of the machine is thus ({I}, {C, K, M }) | I | α 0 . In the following, if the top of a register is #, the register is said empty. We suppose that M is empty at the beginning.

α0 case 0, C jump αC=0 case 1, C jump αC=1 jump αe αe refl M exec M αC=0 pop C case 0, K jump αM=0 case 1, K jump αM=1 jump αe αC=1 pop C case 0, K jump αM=1 case 1, C jump αM=0 jump αe αM=0 pop K move M, 0 jump α0 αM=1 pop K move M, 1 jump α0 Fig. 2. Self-modifying program in ASM2
The program begins by inspecting the first element of the register C (the encrypted message) and decrypts it using the first element in the register K (the key) by jumping to either address α C=0 or address α C=1 . After having been decrypted, it writes the result in the register M (addresses α M=0 and α M=1 ) end goes back to the beginning. If the register C is empty, is goes to the evaluation procedure at address α e which starts the execution of register M .

The self-modifying behavior is clearly exposed thanks to explicit instructions refl and exec.

Further works

Giving this abstract machine, we have first adapted a measure of self-modifying behavior given in [10], the waves, which relates executed codes with codes which wrote them. It appears that this framework is more appropriate to deal with the notion of waves thanks to the time separation of actions of writing and executing, with instruction refl. It makes it possible to more precisely define the moment when a set of addresses can be seen as a program. Annotations make it easier to sequentialize self-modification.

On the other hand, this machine adds a new granularity on the top of addresses, the registers. This allows to structurally see the organisation of the program and infer more easily the intentional behavior of the program. Gathering instructions into high level semantics functions is an important step in the comprehension of programs.

Another field of investigation is to reconstruct self-modifying programs without self-modifications by "unfolding" successive waves of code. This study is based on recovering the executed code of program from execution traces, for any trace. Working with our machine is pleasant because, thanks to separation between data and code, we can more easily define what a program is at a given moment than with x86 semantics where code and data are not distinguishable.

By using this framework, we are also developing a static representation of self-modifying programs based on waves and which takes advantage of larger granularity (registers) of our machine. Indeed, this gives the ability to gather instructions not executed at the same time (for instance with different inputs) hoping that each register contains one and only one high level function. Giving that, we can decide to see statically each program as the tree of all registers which was written and then executed for all inputs, related according to writing actions. This static point of view can be seen as a specification of self-modifying programs, because it describes how each registers are created, that is to say, which register creates another. We can imagine using this specification to design a compilation function taking a non self-modifying program and returning a equivalent selfmodifying program in ASM 2 with respect to the given specification.

  Fig. 1. ASM2 instruction set

⋆ This work is supported by French Direction Générale de l'Armement.