Monika Korte

Richard Holme

On the Persistence of Geomagnetic Flux Lobes in Global Holocene Field Models

Keywords: Geomagnetism, flux lobes, archaeomagnetic field, time-averaged field, millenial secular variation

To investigate the longevity and robustness of high-latitude flux patches in the geomagnetic field at the core-mantle boundary, we present time-dependent models of the geomagnetic field for the past 7000 years. Our models use the same data set as previously used for time-dependent archaeomagnetic field modelling, but constrained with additional priori models from time averages of field models covering the last 150, 400 and 3000 years. We find that the data are consistent with flux patches existing in both north and south hemispheres for the past 7000 years, and that the northern hemisphere patches at least have highly dynamic behaviour. Simple averaging of the historical field may not be appropriate to obtain a characteristic time-averaged model of the field for comparison with other geophysical observables. Our results should inform geodynamo studies of thermal core-mantle coupling, and of possible long-term structure in the geomagnetic field.

Introduction

Near-surface observations of the geomagnetic field provide a powerful probe of the dynamics of the top of the Earth's core, and ultimately of the whole geodynamo. Models of the core surface field have been constructed on a wide range of time scales, from recent, high-resolution models from satellite data (e.g. [START_REF] Lesur | GRIMM: the GFZ Reference Internal Magnetic Model based on vector satellite and observatory data[END_REF][START_REF] Olsen | CHAOS-2 -a geomagnetic field model derived from one decade of continuous satellite data[END_REF], through time-dependent models of the historical [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF] and archaeomagnetic A c c e p t e d M a n u s c r i p t field (e.g. Korte and Constable, 2005;[START_REF] Korte | Geomagnetic field for 0-3ka: 2. a new series of time-varying global models[END_REF], to time-averaged models of the last 5 Myr paleomagnetic field (e.g., [START_REF] Johnson | The time-averaged geomagnetic field as recorded by lava flows over the past 5 Myr[END_REF][START_REF] Kelly | The geomagnetic field over the past 5 Myr[END_REF]. As the time periods increase, data quality and distribution decrease, leading to models with lower resolution in both space and time. Nevertheless, there is evidence of coherent structure in the field on all time scales; such structure is highly significant, and has been interpreted in terms of mantle control on the geodynamo, through thermal [START_REF] Bloxham | Thermal core-mantle interactions[END_REF] or possibly electromagnetic [START_REF] Holme | Electromagnetic core-mantle coupling III: laterally varying mantle conductance[END_REF] core-mantle coupling. Two features of the field have been of particular interest: low secular variation under the Pacific hemisphere, and stationary flux lobes at high latitudes, possibly in symmetric locations in the northern and southern hemisphere. It is on the latter patterns that this paper focuses.

Figure 1 shows the time-averaged field structure at the core-mantle boundary (CMB) on four time scales: the historical period from 1840 -1990 (for which detailed dedicated observations are available, in particular including absolute magnetic intensity determinations, and time-evolution of the field from magnetic observatories), the longer historical period from 1590 -1990 (prior to 1840 dominated by data from ship's logs [START_REF] Jonkers | Four centuries of geomagnetic data from historical records[END_REF]), the archaeomagnetic field for the past 3000 years (inclucing many sedimentary records), and for the past 7000 years. Comparison of the averages shows the expected reduction in resolution with averaging time, particularly in the southern hemisphere. All the averaged global models predict the existence of two or three lobes of strongest magnetic flux in the Northern hemisphere [START_REF] Bloxham | Geomagnetic secular variation[END_REF], with similar features observed in longer term time-averaged global models based on paleomagnetic data from the past 5 million years [START_REF] Johnson | The time-averaged geomagnetic field as recorded by lava flows over the past 5 Myr[END_REF][START_REF] Kelly | The geomagnetic field over the past 5 Myr[END_REF]. However, the southern hemisphere flux patches, clear in recent historical and satellite models, are not seen in the older models.

It is difficult to decide how much field structure we can expect to resolve with archaeoand paleomagnetic data, where data and dating uncertainties are high and often not wellunderstood. Simply truncating the spherical harmonic expansion avoids any small-scale structure, but may map higher degree energy into the lower degree model coefficients. Instead, we seek regularised models: models that both fit the data, and also minimise some [START_REF] Korte | Geomagnetic field for 0-3ka: 2. a new series of time-varying global models[END_REF] and CALS7K.2 (Korte and Constable, 2005), respectively. Northern hemisphere flux patches are seen at high latitudes under North America and Siberia, and in the more recent models approximately symmetric southern counterparts are seen south of South America and Australia.

A c c e p t e d M a n u s c r i p t physical quantity, such as a bound on the electrical dissipation in the dynamo [START_REF] Gubbins | Can the Earth's magnetic field be sustained by core oscillations?[END_REF], the mean square field strength over the CMB, or its time variability. Such assumptions bias the solution for the magnetic field towards the minimum magnitude capable of explaining the observations, and therefore provide a lower bound for the true field strength or complexity. Formally (in a Bayesian sense) we are seeking a model which fits the data, subject to an a priori model of zero magnitude -that there is no field! This assumption is clearly not reasonable: the Holocene data clearly indicate the existence of a persistent field. More specifically, if a time-dependent field model of the archaeomagnetic field suggests episodic flux patches in the north, and none at all in the south, is this because the data demonstrate that the flux patches are not there at some times, or because the model is biased towards low magnitudes, and with limited data, the flux patches are eliminated to reduce the field strength, even (possibly) at the cost of worsened fit to data?

To investigate this issue, we seek models of the field which minimise the deviation of the time-dependent model from a time-averaged model of the field from more recent times, in which the flux patches are seen clearly. Philosphically, an a priori model of a field like that of the present day seems no more unreasonable than a zero prior, which we know to be false. By doing this, we investigate whether it is the data or the regularising assumptions which lead to episodic or missing flux patches. If the models we generate include the flux patches, then we can argue that there is insufficient evidence to reject the patches being consistent features in the field; if the patches move around, are episodic or absent, then we may reject this hypothesis, and make stronger statements about the behaviour of the centenial and millenial structure of the geomagnetic field.

Data and a priori models

The data set is the one used to generate CALS7K.2 (Korte and Constable, 2005), which is readily available from the EarthRef Digital Archive (http://www.earthref.org) (Korte et al., 2005). It consists of directional data from archaeomagnetic studies and lake sediment records, and intensity data determined from archaeological artefacts and lavas. We realize that additional data have been published meanwhile, but in order to allow for a comparison A c c e p t e d M a n u s c r i p t to the CALS7K.2 model we retained that data compilation for this study. An iterative data rejection scheme was used in devising CALS7K.2, and we consider only the final data set of 27067 values. All data which could not be fit within two average standard deviations by a first model had been rejected in this dataset (see Korte and Constable, 2005).

The choice of the prior model is not straightforward. It is not clear over how long a time interval the field has to be averaged in order to minimise rapidly varing smallscale structure to represent properly only persistent larger-scale features. Average over too short a time period, and rapidly varying features are mapped into the stationary mean model; over too long a period and we lose structure due to inadequate data distribution or quality. In order not to make unnecessary assumptions about amount and position of persistent structure we consider three a priori models based on data and without additional filtering: gufm1 [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF] averaged over the time intervals 1840 to 1990 and 1590 to 1990, and the time-average of the recently published 3kyr model CALS3k.3 [START_REF] Korte | Geomagnetic field for 0-3ka: 2. a new series of time-varying global models[END_REF]. All of these averaged models clearly show two pairs of flux lobes, which are approximately hemispherically symmetric. As testing the compatibility of southern hemisphere flux lobes with the 7kyr data set is one of our main motivations, we did not include any longer-term averaged models, which do not show a similarly clear pattern. The motivation for considering gufm1 time averages over both 150 and the full 400 years is twofold: First, 150 years is of the order of the temporal resolution of CALS7K.2 [START_REF] Korte | Spatial and temporal resolution of millennial scale geomagnetic field models[END_REF]. Second, the spatial resolution of this model clearly increases with time as the number and quality of available data increases. In particular, 1840 dates the establishment of geomagnetic observatories, and also start of widespread measurement of absolute intensity. The 150 year average is consequently of higher resolution, containing more smaller scale structure; however, some of this may be present because of insufficient averaging (Fig. 1a andb). The intensity of the model prior to 1840 is unconstrained by data (defined only by a backwards-extrapolation of the dipole strength, cite[]jackson00); this is assumed to be sufficiently close to the real behaviour not to overly bias the average model structure. The gufm1 model is calculated to spherical harmonic degree and order A c c e p t e d M a n u s c r i p t 14; we simply truncated the time-average at degree and order 10 to match the expansion limit of our new models, which seems justified by the fact that with reasonable regularization all the new models show less power than their priors in spherical harmonic degrees of eight and higher.

The CALS3k.3 model is based on an updated archaeomagnetic and sediment dataset [START_REF] Donadini | Geomagnetic field for 0-3ka: 1. new data sets for global modeling[END_REF] and spans the time 1000 BC to 1990 AD. Both spatial and temporal resolution are somewhat higher than for CALS7K.2. The time average (Fig. 1c) contains significantly less structure than the historical averages, but does show clear flux lobes in both the northern and southern hemispheres similar to those seen in the gufm1 averages.

Four new models

The modelling method is basically the same as used for the gufm1 [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF] and CALS7K.2 (Korte and Constable, 2005) models and is described in detail there. The basis functions for the inversion are spherical harmonic functions in space with individual coefficients expanded in cubic B-splines in time to provide a continuous description [START_REF] Bloxham | Time-dependent mapping of the magnetic field at the core-mantle boundary[END_REF]. Maximum spatial and temporal resolution provided by the basis functions are higher than the actual resolution considered feasible from the data with their uncertainties and inhomogeneous global distribution. Regularisations in both time and space are applied in order to find models with minimum structure required by the data.

Many different models may be appropriate solutions to the inverse problem, ranging from very smooth models with large misfit to the data to complex models fitting the data closely.

The regularisation parameters (the damping parameters or Lagrange multipliers controlling the relative penalty assigned to data misfit and model complexity) for the preferred model are commonly found either from the "knee" of a curve trading off misfit against roughness, or by fitting the data to the tolerance given by the uncertainty estimates of the data. For long-term magnetic field models comparisons of resulting time-averaged main field and secular variation geomagnetic power spectra to models constructed with historical data also seem suitable criteria, as higher average spatial or temporal complexity, i.e. higher power in spherical degrees with the exception of 1 and perhaps 2 and 3, is extremely 6 A c c e p t e d M a n u s c r i p t unlikely and certainly not resolvable with the presently available amount and quality of Holocene data [START_REF] Korte | Spatial and temporal resolution of millennial scale geomagnetic field models[END_REF]. The chosen regularistion norms result in a damping of power in main field and secular variation that increases for higher SH degrees, so that small scale / short term structure are efficiently suppressed if the applied regularisation factors result in comparable spectral values for the low degrees [START_REF] Korte | Geomagnetic field for 0-3ka: 2. a new series of time-varying global models[END_REF].

In this study, however, we do not seek an absolute minimum structure model, but the minimum deviation from a given average field structure required by the data. We replace the constraint of minimising a bound on Ohmic dissipation necessary in the dynamo to generate the observed field [START_REF] Gubbins | Can the Earth's magnetic field be sustained by core oscillations?[END_REF], which was used as spatial regularisation for both gufm1 and CALS7K.2, with minimising the radial field deviation at the CMB from an a priori model, i.e. the quantity

S(m 0 ) = CM B (B r -B r (m 0 )) 2 dΩ, (1)
where B r is the radial field of the new model and B r (m 0 ) that of the a priori model m 0 .

The integration is performed over solid angle dΩ at the core surface, averaged over the time period of the model. This condition is easily expressed as a quadratic norm of the geomagnetic Gauss coefficients; we minimise

lmax l=1 (l + 1) 2 2l + 1 a c 2l+4 l m=1 (g m l -g m l (m 0 )) 2 + (h m l -h m l (m 0 )) 2 (2)
where a is the radius of the Earth, c the radius of the CMB, {g m l , h m l } are the geomagnetic Gauss coefficients of spherical harmonic degree l and order m, and {g m l (m 0 ), h m l (m 0 )} the coefficients of the a priori model.

Studies of virtual axial dipole moment [START_REF] Mcelhinny | Variations in the geomagnetic dipole: I. The past 50 000 years[END_REF][START_REF] Yang | Variations in the geomagnetic dipole moment over the last 12000 years[END_REF] and the previous millennial scale models indicate that the dipole moment has varied significantly over the past 7 kyrs. The dipole is the strongest field contribution; we were concerned that taking it into account in investigating the required deviation might have a dominating influence. Therefore, we also tested models where either the axial dipole coefficient, or all three dipole coefficients were not influenced by the spatial regularisation.

A c c e p t e d M a n u s c r i p t

We found that with our criteria for the preferred amount of regularization the differences between these three types of model were small. Nevertheless we retained only the models where the dipole was not included in the spatial regularization in the following comparison.

The final modelling procedure minimises the functional

RMS 2 + λ S S(m 0 ) + λ T T (3)
with spatial and temporal damping factors (Lagrange multipliers) λ S and λ T respectively.

The normalised root mean square misfit (RMS) between model predictions xi and data x i is defined as

RMS = 1 N N i=1 x i -xi σ i 2 , (4)
with uncertainty estimates σ i and N the number of data. The spatial norm compared with the a priori model m 0 is defined in equation ( 1), and the temporal norm (like the spatial norm, averaged over the modelling period) is defined

T = ∂ 2 B r ∂t 2 2 dΩ. (5)
(t is time) which can be calculated using a quadratic norm of the form of equation ( 2).

In the end we calculated four new models, summarized in tabletab:models.

Results

As a first step to comparing our models, Figure 2 shows the main field and secular variation spectra of the time averages of the four investigated models at the CMB together with the spectra of the a priori models and the CALS7K.2 spectra. By plotting these spectra at the CMB we emphasize the differences in higher degrees. The secular variation spectra of the four new models look nearly identical, which is not surprising due to the common temporal damping applied independently from the spatial damping.

The time-averaged spectra of the new models clearly resemble those of their respective a priori models, except in the case of M0 which was designed to be similar to CALS7K.2.

The increasing deviation with spherical harmonic degree above degree l = 5 between M0 and CALS7K.2 is due to the different regulariation applied.

We next consider global diagnostics averaged over the whole 7 kyr time interval. We calculate the RMS misfit (equation 4), departure from a priori model S(m 0 ) (equation 1)

and temporal norm T (equation 5). We also consider the overall spatial structure S(0)

S(0) = CM B B 2 r dΩ. (6)
For M0, with the zero a priori model, we obviously have S(m 0 ) = S(0). These diagnostics are given in Table 2. (Note that although the dipole terms are excluded from the damping, all of the norm calculations include the dipole components to give quantities that have di- In conclusion, the data can be fit better by a model requiring less temporal variability if a suitable a priori model is used. Among the a priori models tested the averaged CALS3k.3 turns out to be the most suitable, with relatively little deviation required, yet giving the best fit to the data and least required temporal variation. The data are least compatible with gufm1 averaged only over the most recent 150 years.

A more detailed analysis of the behaviour of the solution norms as a function of time is provided in Fig. 3, which plots the mean square radial field, the mean square departure of the field from the a priori model, and the temporal norm as a function of time. The profiles for the different a priori models are generally similar, especially for the temporal complexity T . The absolute amount of complexity, S(0), shows nearly identical relative variations with time. The relative deviation from the three non-zero a priori models also is similar, with maximum values around 4000 BC for all, and again between 0 and 1000 AD for the two gufm1 a priori models.

The average spatial structure of the models is shown by plots of their mean radial field at the CMB (Fig. 4). The averages of the different models show clear similarities in large-scale features, while reflecting the amount of complexity of the a priori model (for example, the small near-equatorial flux patches in M150 are clearly a result of insufficient time-averaging of the a priori model in this region). The available data are clearly compatible with two southern hemisphere flux lobes which are persistent enough to show up in the time-averaged model. In all four models, there are three (rather than two) flux lobes present in the northern hemisphere, despite the prior model for M150 and M400 requiring only the two lobes seen in today's field. We may therefore be confident that the Holocene data do require a third northern flux lobe under Europe. A similar flux lobe pattern could also exist in the southern hemisphere, but none of the prior models show this feature, and there are insufficient data to constrain this question. Note, however, that even

14

A c c e p t e d M a n u s c r i p t hemisphere). This apparent third flux lobe therefore might be a manifestation of unresolved but non-averaging smaller scale structure. Another common feature of all models is an area of positive radial field in the north-western Pacific, although its detailed form and strength is affected by the different prior models. This feature has also been seen in previous models (see, for example, the snapshot models of [START_REF] Constable | Global geomagnetic field models for the past 3000 years: transient or permanent flux lobes? Phil[END_REF]), and shows up clearly in longer time-period palaeomagnetic models (for example, [START_REF] Johnson | The time-averaged geomagnetic field: global and regional biases for 0-5 Ma[END_REF].

Animations of the evolution radial field are provided in

Discussion

Our primary result is clear and perhaps unsurprising: southern hemisphere flux patches are consistent with the available data. There is no evidence requiring that they are less persistent then their northern hemisphere counterparts. However, some interesting further results emerge from more detailed comparisons. Constraining the model about the 150 year time average is apparantly less appropriate than applying no constraint at all. This suggests that, although tempting because of the much higher quality data for this period, using this time average as a proxy for long-term field behaviour is not appropriate. Why might this be? One possibility is that the averaging time is insufficient to average out small scale motions (for example, the propagation of flux patches along the equator [START_REF] Jackson | Intense equatorial flux spots on the surface of the Earth's core[END_REF]), leaving small scale features to be fit that are not persistent on longer time scales. Another important effect is likely to be the appearance after 1840 of the southern hemisphere reversed flux patch, associated by [START_REF] Gubbins | Mechanism for geomagnetic polarity reversals[END_REF] with the current rapid decay of the axial dipole field. A third effect is that the particularly strong flux lobes from the a priori model for M150 act to preclude the variability in the position of the northern hemisphere flux lobes; from examination of the movie, there is certainly evidence that the prior model is having a strong influence on the position and longevity of these features compared with the other three models. Any or all of these explanations suggest that if we wish to correlate geomagnetic field morphology with other geophysical observables (e.g. geoid, seismic tomography), then the recent historical field or even high-resolution models of the current field from satellite data may be less appropriate than the longer term averages.

It is interesting to note that all three models using a priori field averages for regular- mainly show two (Fig. 1). The highly dynamic evolution of the flux pattern over time, however, makes it difficult to clearly distinguish between lateral movements of flux patches and growth and decay of regional flux concentrations. Interestingly, a strong appearance of all three flux patches simultaneously is rare. Their movement with time is also significant, because the prior models (M150 in particular) favour restraining them to a single location;

that they are nonetheless variable in their positions suggests that this is a true feature of the field required by the data. Nevertheless, the dynamic nature of the flux patches has important implications for numerical dynamo studies of thermal core-mantle example.

For example, in their dynamo calculations, [START_REF] Willis | Thermal core-mantle interaction: Exploring regimes for 'locked' dynamo action[END_REF] have located a region of parameter space (admittedly far from Earth-like) in which a numerical dynamo code yields flux lobe patterns similar to field observed for the present day. These patches are dynamic, moving around and occasionally dividing, but none the less, apparantly less dynamic than the behaviour implied by the observations and our modelling here. Additional study [START_REF] Davies | Time-averaged paleomagnetic field and secular variation: Predictions from dynamo solutions based on lower mantle seismic tomography[END_REF] has located a parameter regime with evidence of the three-fold symmetry we observe in our models; our models suggest that an evaluation of the transition between these two regimes would be of great interest.

No corresponding southern hemisphere counterpart for the third flux lobe centered under Europe appears. This is not surprising as our a priori models do not encourage this, and the southern hemisphere data are sparse. However, the three models do suggest a region of weak to reversed flux under southern Africa and surrounding areas, roughly the same region as the present day Southern Atlantic Anomaly [START_REF] Gubbins | Geomagnetic field ananlysis -III. Magnetic fields on the core-mantle boundary[END_REF]. Together with the appearance of a strong reverse flux patch around Southern Africa in the first two millennia of the models, this might be interpreted as a preferred area for the recurrence of significant reverse flux and consequently minimum field strength. However, this southern hemisphere feature results from the strong inclination variation seen in one African sediment record between 4000 and 5000 BC (Lake Victoria by [START_REF] Mothersill | Paleomagnetic results from lakes Victoria and Albert, Uganda[END_REF]).

The reverse flux area south of the equator seen in the averages in Fig. 4 disappears if this specific sediment record is omitted from the modelling. To our knowledge there are to

Conclusions

We have investigated persistent structure in the time-averaged geomagnetic field on time-scales from centuries to millennia. The deviation from different time-averaged field models has been used as regularization constraint for spatial structure in modelling the 7kyr dataset which had previously been used for the CALS7K.2 model. Comparisons of misfit and temporal variability resulting with fixed modelling parameters show that the data are less compatible with field averages of the past 150 or 400 years than with a 3kyr average. This indicates that small-scale structure present in field averages of a few centuries is not persistent on longer time-scales. The smallest misfit between data and model is obtained, however, if a 3kyr average instead of a zero assumption is used as smoothing constraint. Distinct northern and southern hemisphere flux lobes are clearly compatible with the available data spanning the past 7kyrs. Note, however, that persistence in time averages of the field does not exclude significant temporal variations on shorter intervals; indeed, from the detailed temporal behaviour there is evidence for considerable variability in these features.

However, clever modelling can take us only so far. While we have demonstrated that currently available data do not preclude southern hemisphere flux patches, only by expanding the data base for this region can we truly determine the long-term morphology and variability of the southern hemisphere magnetic field.

Figure 1 :

 1 Figure 1: Time averaged radial field component at the core-mantle boundary of the three a priori models and CALS7K.2. (a) and (b) are 150 and 400 year averages, respectively, of gufm1 by Jackson et al. (2000), (c) and (d) are the overall time averages of CALS3k.3 (Korte et al., 2009) and CALS7K.2 (Korte and

  For comparability, we first reconstructed a model M0 similar to CALS7K.2 with zero prior model, but regularising by minimising the mean square radial field at the CMB instead of Ohmic dissipation norm, and also excluding the dipole coefficients from the regularization. We chose the damping factors such that the resulting model shows similar main field and secular variation spectra to CALS7K.2 and has a comparable data misfit -the values are λ S = 5 × 10 -11 and λ T = 10 2 . Three further models were constructed with the same damping parameters; the only difference being the a priori model. The averaged 150 and 400 year gufm1 were used respectively for models M150 and M400. CALS3k.

Figure 2 :

 2 Figure 2: Time averaged main field (top) and secular variation (bottom) geomagnetic power spectra at the core-mantle boundary for models M0 (open triangles), M150 (open squares), M400 (open circles), M3k (open diamonds) and time averaged gufm1 for 1840 to 1990 (gray squares), 1590 to 1990 (gray circles), CALS3k.3 (gray diamonds) and CALS7K.2 for 5000 BC to 1950 AD (gray triangles).

  M400 marginally better, and the M3k substantially better.

  recent high-resolution field models for one epoch, like e.g. the IGRF for 2005 (Macmillan and Maus, 2005), show a somewhat similar third flux lobe if truncated at spherical harmonic degree 5 or 6 (with counterpart in the southern A c c e p t e d M a n u s c r i p t

Figure 3 :

 3 Figure 3: Spatial and temporal complexity as measured by S(m 0 ) (black solid line), S(0) (dashed line) and T (gray) for models a) M0, b) M150, c) M400, d) M3k. The solid and dashed lines are by definition identical for M0.

Figure 4 :

 4 Figure 4: Time averaged radial field component at the core-mantle boundary of the four new models.

Fig. 5 (

 5 electronic version or supplemental material) with the present position of flux lobes outlined by the +/-400 µT isolines of the 400 yr time-averaged gufm1 model. These animations show substantial variability of the flux concentrations on multi-centennial time-scales in all four models, confirming that they are in fact required by the data and with only minor influence from the choice of a priori model. Despite a large amount of movement and/or decrease and increase of flux, the flux concentration is rather high most of the time in the area of the present North American flux lobe, with only two time intervals of significantly weaker flux spanning about 500 years around 1950 BC and 650 AD. Similarly, the flux concentration remains high in the area of the present Siberian flux lobe except for 200 to 600 year intervals around 3450 BC, 1800 BC and 100 BC. Significantly stronger flux than present, however, appears in the European region for 2 to 5 centuries around 5000 BC, 4500 BC, 2500 BC, 1500 BC, 250 BC and 800 AD. Flux variations are weaker in the southern hemisphere, likely as a consequence of the sparcity of southern hemisphere data.

Figure 5 :

 5 Figure 5: Animations (electronic version or supplemental material) of the radial magnetic field evolution of models (a) M0, (b) M150, (c) M400 and (d) M3k at the CMB. Yellow lines indicate the +/-400 µT contour lines of the 400 yr time averaged gufm1 model.

  ization cleary indicate the existence of three large northern hemisphere areas of high flux Page 17 of 25 A c c e p t e d M a n u s c r i p t concentration on average (Fig. 4), while the historical field averages (and present field)

  to the south of that location between Argentina and Australia to support or contradict the strong effect of these data on the model. Until such data are available this model feature should be regarded with caution.

Table 1 :

 1 The four new models and their priors.

	Model Prior
	M0	0
	M150 gumf1 averaged from 1840 to 1990
	M400 gufm1 averaged from 1590 to 1990
	M3k	CALS3k.3 averaged over 3kyr

Table 2 :

 2 Rms misfit and spatial and temporal complexity of new modelsModelrms T (nT 2 yr -4 ) S(0) (nT 2 ) S(m 0 ) (nT 2 )

	Time interval 5000 BC to 1950 AD:
	M0	0.989	10.69 5.9 × 10 10	5.9 × 10 10
	M150	0.997	10.85 7.6 × 10 10	2.0 × 10 10
	M400	0.990	10.63 7.0 × 10 10	1.6 × 10 10
	M3k	0.984	10.53 6.3 × 10 10	1.2 × 10 10
	Time interval 5000 BC to 1000 BC:
	M0	0.968	7.87 4.8 × 10 10	4.8 × 10 10
	M150	0.975	7.94 6.6 × 10 10	3.0 × 10 10
	M400	0.961	7.62 6.0 × 10 10	2.3 × 10 10
	M3k	0.956	7.49 5.3 × 10 10	1.6 × 10 10

dataset for the time interval 1000 BC to 1950 AD form a subset of the data used in the construction of CALS3k.3, the a priori model for M3k. This will inevitably give a bias in our results; the time-averaged CALS3K model is bound to be a good fit to the CALS7K data set for the period after 1000 BC. To accommodate this bias, we calculate two sets of diagnostics in table 2, first for the whole 7000 years, and secondly separately for the time interval 5000 BC to 1000 BC, which is not considered in any of the a priori models. This second set of diagnostics should eliminate (or at least limit) the bias from considering common data sets. With so many diagnostics of the models presented, it could be extremely difficult to determine which prior model is most compatible with the observations; expected). Comparing misfit and temporal norm, the M150 model performs less well than

Acknowledgements

We thank two anonymous reviewers who provided very constructive comments that improved the mansucript. Funding for RH and the collaboration which produced this research was provided by NERC grant NER/O/S/2003/00675. The project was completed with support of a DAAD travel grant.

A c c e p t e d M a n u s c r i p t


J Nisar

R Ahuja

Equation of state (EOS) and collapse of magnetism in Iron-rich meteorites at high pressure by first-principles calculations

Keywords: Iron-rich composition, Equation of state, Collapse of magnetism

We have calculated electronic and structure properties, equation of state (EOS) and collapse of magnetism of different meteoritic materials at high pressure using density functional theory (DFT). The fully optimization structures are obtained by minimization of the total energy and atomic forces of iron rich materials, which are in good agreement with the experiment. Phase 1.

transition of Fe 2 P (C22 type structure) barringerite into Fe 2 P (C23 type structure) allabogdanite is calculated to be 26 GPa. Collapse of magnetism (high spin-low spin) of some iron rich materials are described at various pressures. It has been observed that the optimized equilibrium volume has been changed due to transition of non-magnetic to magnetic material at high pressure, so EOS parameters have also been changed. EOS parameters such as bulk modulus (B o ) and derivative of bulk modulus with respect to pressure ( o B′ ) for magnetic and nonmagnetic configurations have been calculated by fitting E-V curves with Birch-Murnaghan equation, which confirm the experimental results. Introduction

High pressure study of materials has great importance for the understanding of phase transition mechanism, mechanical stability and the internal structure of Earth core and other planets. The Bulk modulus of material determines mechanical properties of the material and the amount of energy stored in solid material in the Earth's core. The equation of state (EOS) is an important parameter for characteristic of a material. One of the most abundant constituent of the Earth core is iron. The Earth's outer core also contains certain amount of light elements because the difference in the density between solid iron and the inner core may be as much as 5%.

( [START_REF] Mcdonough | The composition of the Earth[END_REF][START_REF] Anderson | An equation of state for liquid iron and implications for the Earth's core[END_REF][START_REF] Sohl | The interior structure of Mars: implication from SNC meteorites[END_REF]. Phosphorous, carbon and sulpur are believed to be alloying light elements in iron-rich planetary cores, such as, those of the Earth and Mars [START_REF] Li | Experimental constrains on core composition[END_REF][START_REF] Poirier | Light elements in the Earth's core: a critical review[END_REF].

Determination of the physical properties of the magnetic and non-magnetic phase is important for the interpretation of existence of light elements in Earth and planetary cores. The study of density difference between pure iron and light element contained compounds such as Fe 2 P, Fe 3 P, Fe 3 S and Fe 3 C is very important for the estimation of light elements in Earth's core.

Physical properties of meteoritic materials under core conditions can be significantly affected by the state of magnetization. High pressure magnetism in transition metal oxides had been investigated with in density function theory (DFT) [START_REF] Cohen | Magnetic collapse in transition metal oxides at high pressure: Imlications for the Earth[END_REF]. Fe 2 P is a very common material which is found in two structures as meteoritic samples. Hexagonal Fe 2 P barringerite with space group m P 2 6 _ (C22 structure type) is found in ambient pressure and orthorhombic polymorph of Fe 2 P allabogdanite with space group Pnma (C23 type structure) are also known as stable or metastable at ambient conditions [START_REF] Britvin | Allabogdanite, (Fe, Ni) 2 P, a new mineral from the Onello meteorite: The occurrence and crystal structure[END_REF].

A c c e p t e d M a n u s c r i p t

In this present work, the lattice parameters of expected Earth core materials are measured. We have used first principle calculation to determine EOS parameters and volume in both the magnetically ordered and non-magnetic configurations of Fe 2 P (barringerite), Fe 2 P (allabogdanite), Fe 3 P, Fe 3 S, Fe 3 C, Fe (bcc), and Fe (hcp). The magnetic transition (high spin to low spin) pressures of iron rich materials have been calculated. Magnetic transition in the material can occur due to the widening of bands [START_REF] Ohnishi | A theory of the pressure-induced high spin -low spin transition of transition metals oxides[END_REF]. Phase transition of Fe 2 P (barringerite) in Fe 2 P (allabogdanite) has also been investigated in this paper.

Computational Details

We have studied the high-pressure-induced electronic structure, EOS, magnetic transition in iron-rich materials by performing first-principles calculations based on the DFT [START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF]. Ab initio simulations were performed using the projector-augmented wave (PAW) method [START_REF] Blochl | Projector augmented-wave method[END_REF] as implemented in Vienna ab initio simulation package (VASP) code [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmentedwave method[END_REF]. Generalized gradient approximation (GGA) (Perdew et al., 1992) was used to describe the electronic exchange-correlation effects. Fe 3 P and Fe 3 S have tetragonal structure with space group _ 4 I , and Fe 3 C has orthorhombic structure with space group Pnma. We used denser K-point grid for total energy calculations. The calculations were carried out using 7×7×11, 9×9×9, 8×6×8, 16×16×16, 9×9×7, 3×3×7 and 3×3×7 k-points meshes for Fe 2 P (C22), Fe 2 P (C23), Fe 3 C, Fe(bcc), Fe(hcp), Fe 3 P and Fe 3 S, respectively. All Brillouin zone integrations were performed with a Monkhorst-Pack generated k-points mesh [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF]. The PAW potentials have been used with 4s, 3d electrons in valence for iron, 3s, 3p electrons for phosphorous and sulfur and 2s, 2p electrons in valence for carbon. A kinetic energy cutoff of 400 eV was used for the plane wave basis set for all materials. The atomic geometries were fully optimized until the Hellmann-Feynman forces on each atom were less than the threshold value of 10 -4 eV/Å. For the A c c e p t e d M a n u s c r i p t Brillouin zone integrations, the Gaussian smearing method with a 0.02 eV smearing width has been used for the geometry optimization.

In DFT, the interaction between valence and core electrons are described by pseudo potentials and valence orbital which are expressed as plane waves. All calculations have been performed at 0 K and vibrational contributions are not included. The EOS is calculated by using the internal energy of the crystal at different chosen volume. The resulting energy versus volume data have been fitted to an integrated third-order Birch-Murnaghan (BM) equation [START_REF] Birch | Finite strain isotherm and velocities for single-crystal and polycrystalline sodium chloride at high pressures and 300[END_REF], which is given in equation 1.

( ) ( ) ( ) ( ) ( ) 1 3 16 3 14 4 16 9 3 2 3 5 3 4 3 7 2 3 o o o o o o o o E V V B V V B V V B B V E + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ′ - + ′ - - ′ - - =
The pressure is obtained by taking volume derivative of the total energy. The bulk modulus B o = -V o dP/dV is also calculated from P-V relation where V o is the equilibrium volume obtained from full structural relaxation and o B′ is the first derivative of bulk modulus with respect of pressure.

Both magnetic (spin polarized with ferromagnetic ordering) and non-magnetic (non-spin polarized) treatment have been considered. The total energy has been obtained by fully selfconsistent ab initio electronic structure calculations for each configuration.

Results and Discussion

In table 1, the structure parameters of some meteoritic materials and the equilibrium unit cell volume for each structure obtained by geometrical optimization procedure using firstprinciples method, which are very close to the experimental results. Deviation between experimental and theoretical values of unit cell volume is varied from 2 to 3.6 % (Fe 2 P (C22), Fe 2 P (C23), Fe 3 P, Fe 3 C, Fe 3 S, Fe (bcc) and Fe (hcp)).

A c c e p t e d M a n u s c r i p t

Sulpur, phosphorous and silicon are the main candidates for the light element core constituent. Fe 2 P is ferromagnetic in both crystal structure and structural phase transformation, which occur from hexagonal barringerite with the C22 type structure to orthorhombic allabogdanite with the C23 type structure under high pressure (Senateur et al., 1976). Our calculations describe that hexagonal Fe 2 P undergoes structure phase transformation into orthorhombic Fe 2 P at 26 GPa. We have compared the total density of states (TDOS) of orthorhombic Fe 2 P and hexagonal Fe 2 P in Fig. 1, which also indicates that Fe 2 P (orthorhombic)

is more stable at 26 Gpa. Both structures are ferromagnetically stable but magnetic moment per Fe atom of Fe 2 P (C22) is higher than Fe 2 P (C23). The calculated phase transition pressure of Fe 2 P is somewhat higher than the experimental value (Dera et al., 2008). Main reason of this theoretical overestimation is the temperature, because experimental phase transformation is observed at 8 GPa and 1400 K.

Both magnetic and non-magnetic calculations were performed for some iron rich alloys.

The total energy and the volume curves of different materials were fitted by third order BM equation [START_REF] Anderson | Equations of state of solids for Geophysics and Ceramics Science[END_REF] to obtain the bulk modulus (B o ), its first derivative with respect to pressure ( o B′ ) and equilibrium volume for each structure. The results are presented in Table 2 together with the results from other theoretical prediction and available experimental data. The non-magnetic phase of materials has higher incompressibility and density than the magnetic phase. We find that there is a significant difference in the EOS parameters of the magnetic and nonmagnetic phase of materials. Studies of the non-magnetic phase of iron rich materials are more relevant than the magnetic phase for understanding the properties of materials under the conditions of planetary cores because of high pressure and temperature. Bulk modulus (B o ) and its first derivative ( o B′ ) of Fe 3 C for non-magnetic calculations are 297.3 GPa and 4.9 respectively 2002). Magnetic (spin polarized) calculation of EOS parameter such as bulk modulus of some meteoritic materials are slightly overestimated as compared with the experimental results.

The magnetic collapse (high spin-low spin) of materials have been found due to the contribution of substantial volume reduction and shortening of ionic distances as described in the study of high pressure magnetic transition of hermatite (Fe 2 O 3 ) [START_REF] Rozenberg | High-pressure structural studies of hematite Fe 2 O 3[END_REF][START_REF] Badro | Nature of the high-pressure transition in Fe 2 O 3[END_REF]. Fig. 2 shows the calculated variation of magnetic moment with pressure of some meteoritic materials. It is clear from the figure that magnetic moment of meteoritic materials decrease with the increase of pressure, which is due to band widening. The magnetic to nonmagnetic phase transition affects the elastic, thermodynamics, compressibility and mechanical properties of materials. We have observed that pressure induced magnetic transition (high spin to low spin) of Fe 3 C and Fe 3 S at pressures of 32-38 and 24-30 GPa respectively, which are slightly larger than the experimental values 25 GPa , 21 GPa respectively (Lin et al., 2004, Lin et al., EPSL 2004). It is well known that magnetic moment of material reduces or become low at higher temperature. Our DFT calculations are performed at 0 K but all experiments are done at 300 K.

This could be one of the reasons that our results are always overestimated. Collapse of magnetism of Fe 2 P (C23) occurs at 56 GPa approximately. Magnetic transitions of iron rich materials affect the description of the amount and properties of light elements in the Earth's core.

Volume differences between magnetic and non-magnetic calculation of iron-rich compounds are between 2.8 to 7.42 % (see in table 2). Physical properties such as compressibility and thermal expansion of the iron rich materials are changed due to volume reduction at magnetic transition pressures.

A c c e p t e d M a n u s c r i p t

In this work, Ferromagnetism of iron rich materials has also been investigated.

Ferromagnetic stability has been determined by the total energy difference (∆E = E AFM -E FM )

per unit cell between the anti ferromagnetic (AFM) configuration and ferromagnetic (FM) configuration. If ∆E is negative, the AFM configuration is more stable than FM and vice versa.

At ambient pressure, materials such as Fe 2 P (C22), Fe 2 P (C23), Fe 3 P, Fe 3 S, Fe 3 C and Fe (bcc) are ferromagnetic stable with average magnetic moment 1.4, 0.7, 1.8, 1.9, 1.7 and 2.2 µ B per iron atom, respectively. Our results agree with previous experimental and theoretical results (Vocadlo et al., 2002[START_REF] Caspersen | Importance of Shear in the bcc-tohcp Transformation in iron[END_REF].

Conclusions

To summarize the complete analysis of the structure, electronic and magnetic properties of some meteoritic materials have been presented using first-principles technique. EOS parameters of some Earth core materials (Fe 2 P (C22), Fe 2 P (C23), Fe 3 P, Fe 3 S, Fe 3 C and Fe (bcc)) are calculated for both magnetic and nonmagnetic calculations, which are in good agreement with experimental results but some of them are slightly overestimated. The firstprinciples calculations also show that the non-magnetic phase has a high bulk modulus than magnetic phase. It is also observed that phase transition of Fe 2 P (barringerite) into Fe 2 P (allabogdanite) is at 26 GPa. Collapse of magnetism of these meteoritic materials is described at different pressures. The observed magnetic transition would affect the interpretation of the amount and properties of the light elements present in the planetary cores. All above mention meteoritic materials are ferromagnetic stable at ambient pressure except Fe (hcp), which is antiferromagnetically stable. We hope our work will prompt more experimental studies about magnetism under high pressure for these materials. c Fei et al., (2000) d [START_REF] Mao | Effect of pressure on crystal structure and lattice parameters of iron up to 300 kbar[END_REF] A c c e p t e d M a n u s c r i p t

  good agreement with previous results, Bo = 316.6 and o B′ =4.3 (Vocadle et al.,

Fig. 1 .

 1 Fig. 1. Comparison of total density of state (TDOS) of Fe 2 P (C22) and Fe 2 P (C23) at 26 GPa, the

Table 1 .

 1 Optimized Lattice parameters (a, b and c, Å) and cell volumes (V, Å 3 ) of iron-rich materials

			This work (GGA)			Experimental
	Materials
		a	b	c	V	a	b	c	V
	Fe 2 P (C22)	5.794	-	3.440	100.01	5.876 a	-	3.449	103.15
	Fe 2 P (C23)	5.543	3.559	6.503	128.32	5.381 a	3.447	6.422	119.12
	Fe 3 P	9.027	-	4.369	355.96	9.107 b	-	4.46	369.9
	Fe 3 S	9.08	-	4.348	359.32	9.144 c	-	4.509	377.0
	Fe 3 C	4.995	6.759	4.432	149.67	5.058 d	6.703	4.506	152.77
	Fe (bcc)	2.83	-	-	22.63	2.866 d	-	-	23.54
	Fe(hcp)	2.445	-	3.886	20.1	2.523 d	-	4.037	22.25

a

Dera et al., (2008)

b

Rundqvist., (1962)

Table 2 .

 2 Equation of state (EOS) parameters for different light element containing components of Earth core (Bulk modulus, Bo, in GPa and its first derivative with respect to pressure, o B′ ) and equilibrium volume, V o , Å 3 )

	Material	Non-spin Polarized Calculation		Spin Polarized Calculation		Previous data
		Bo	o B′	V o	Bo	o B′	V o	Bo

Acknowledgment

We would like to acknowledge VR for financial support. One of us, Jawad Nisar also acknowledges Higher Education Commission (HEC) of Pakistan for financial support. SNIC and UPPMAX are acknowledged for providing computing time.


Dimensionality analysis of magnetotelluric data is a common procedure for inferring the main properties of the geoelectric structures of the subsurface such as the strike direction or the presence of superficial distorting bodies, and enables the most appropriate modeling approach (1D, 2D or 3D) to be determined. Most of the methods currently used assume that the electrical conductivity of individual parts of a structure is isotropic, although some traces of anisotropy in data responses can be recognized. In

this paper we investigate the imprints of anisotropic media responses in dimensionality analysis using rotational invariants of the magnetotelluric tensor. We show results for responses generated from 2D synthetic anisotropic models and for field data that have been interpreted as showing the effects of electrical anisotropy in parts of the subsurface structure. As a result of this study we extend the WAL dimensionality criteria to include extra conditions that allow anisotropic media to be distinguished from 2D isotropic ones. The new conditions require the analysis of the strike directions obtained and take into account the overall behavior of different sites in a survey.

Introduction

Electrical anisotropy in the Earth, caused by electrical conductivity varying with orientation, is a property that is increasingly being taken into account in the interpretation of magnetotelluric data. Electrical anisotropy in the crust can be caused by preferred orientations of fluids, sulfides or fractures [START_REF] Wannamaker | Anisotropy versus heterogeneity in continental solid Earth electromagnetic studies: Fundamental response characteristics and implications for physicochemical state[END_REF], whereas in the upper mantle, it is linked to the splitting of seismic SKS waves [START_REF] Eaton | Tectonic fabric of the subcontinental lithosphere: Evidence from seismic, magnetotelluric and mechanical anisotropy[END_REF], and is explained by either hydrogen diffusivity in olivine crystals [START_REF] Wannamaker | Anisotropy versus heterogeneity in continental solid Earth electromagnetic studies: Fundamental response characteristics and implications for physicochemical state[END_REF][START_REF] Wang | The effect of water on the electrical conductivity of olivine[END_REF] or by the presence of partial melt elongated in the direction of plate motion [START_REF] Yoshino | Hydrous olivine unable to account for conductivity anomaly at the top of the asthenosphere[END_REF].

Significant developments have been achieved regarding the study of electrical anisotropy using magnetotellurics. These deal with modelling and inversion schemes, which include anisotropy [START_REF] Pek | Finite-difference modelling of magnetotelluric fields in twodimensional anisotropic media[END_REF][START_REF] Weidelt | 3D conductivity models: implications of electrical anisotropy[END_REF][START_REF] Wang | 3-D electromagnetic anisotropy modeling using finite differences[END_REF][START_REF] Li | A finite-element algorithm for electromagnetic induction in twodimensional anisotropic conductivity structures[END_REF][START_REF] Yin | Inherent nonuniqueness in magnetotelluric inversion for 1D anisotropic models[END_REF]Pek andSantos, 2002, 2006), the analysis of magnetotelluric responses affected by anisotropy [START_REF] Reddy | Magnetotelluric response of laterally inhomogeneous and anistotropic media[END_REF][START_REF] Saraf | Magnetotelluric response of a laterally inhomogeneous anisotropic inclusion[END_REF][START_REF] Osella | Magnetotelluric response of anisotropic 2-D structures[END_REF][START_REF] Heise | Anomalous phases exceeding 90º in magnetotellurics: anisotropic model studies and a field example[END_REF][START_REF] Heise | Anisotropy and phase splits in magnetotellurics[END_REF], and the investigation of the intrinsic properties and processes causing electric anisotropy [START_REF] Gatzemeier | Flow and electrical anisotropy in the upper mantle: Finite-element models constraints on the effects of olivine crystal preferred orientation and microstructure[END_REF]. Some of the aforementioned papers were published in a special issue dedicated to electrical and seismic continental anisotropy [START_REF] Eaton | Tectonic fabric of the subcontinental lithosphere: Evidence from seismic, magnetotelluric and mechanical anisotropy[END_REF]. A review of earlier work can be found in [START_REF] Wannamaker | Anisotropy versus heterogeneity in continental solid Earth electromagnetic studies: Fundamental response characteristics and implications for physicochemical state[END_REF].

To date there have been no studies specifically discussing the effects of anisotropy on rotational invariants or its complete dimensionality characterization. The goal of this paper is to identify electrical anisotropy using dimensionality analysis based on the rotational invariants of the magnetotelluric tensor. Data were generated from various A c c e p t e d M a n u s c r i p t 4 synthetic models with electrical anisotropy using the 2D code of [START_REF] Pek | Finite-difference modelling of magnetotelluric fields in twodimensional anisotropic media[END_REF].

The results from a set of field data that has been interpreted as exhibiting the effects of anisotropic Earth structure (from the COPROD dataset) are also discussed.

Background

Dimensionality analysis in magnetotellurics

In the magnetotelluric (MT) method (e.g. [START_REF] Vozoff | The magnetotelluric method[END_REF][START_REF] Simpson | Practical Magnetotellurics[END_REF], dimensionality analysis is a common procedure for determining, prior to modeling, whether the measured data or computed responses (impedance tensor, Z; tipper, T; apparent resistivities, ρ ij ; and phases, φ ij ) at a given frequency (ω) correspond to 1D, 2D or 3D geoelectrical structures. It also allows identification and quantification of distortions [START_REF] Kaufman | Reduction of the geological noise in magnetotelluric soundings[END_REF][START_REF] Groom | Decomposition of the magnetotelluric impedance tensor in the presence of local three-dimensional galvanic distortion[END_REF][START_REF] Smith | Understanding telluric distortion matrices[END_REF] and, when applicable, recovery of the directionality (strike) of the structures. Dimensionality analysis techniques search for particular relationships between the components of the magnetotelluric impedance tensor, Z(ω) (e.g. [START_REF] Cantwell | Detection and analysis of low frequency magnetotelluric signals[END_REF], or related functions, in order to identify each dimensionality type. Additional information can be obtained from the induction arrows (i.e., tipper vectors). The dimensionality analysis technique that sees the most widespread use is that of [START_REF] Mcneice | Multisite, multifrequency tensor decomposition of magnetotelluric data[END_REF]. This technique uses the [START_REF] Groom | Decomposition of the magnetotelluric impedance tensor in the presence of local three-dimensional galvanic distortion[END_REF] decomposition method to find the best fitting 2D parameters for a set of sites at different period bands. [START_REF] Lilley | Magnetotelluric analysis using Mohr circles[END_REF] introduced the use of Mohr circles to display and analyze magnetotelluric data, allowing to distinguish their dimensionality and the presence of galvanic distortion. In two-dimensional cases, the regional geologic strike is estimated from either the real or imaginary parts of the magnetotelluric tensor ( hr ,  hq , eqs. 113 and 114, Lilley, 1998a) [START_REF] Weaver | Characterisation of the magnetotelluric tensor in terms of its invariants[END_REF] (based on [START_REF] Lilley | Magnetotelluric analysis using Mohr circles[END_REF]Lilley, , 1998a;;Fisher and Masero, 1994;[START_REF] Szarka | Analysis of rotational invariants of the magnetotelluric impedance tensor[END_REF]) presented a complete dimensionality criteria based on the rotational invariants (WAL invariants) of the magnetotelluric tensor ( () M  , defined as the relationship between the electric field E(ω) and the magnetic induction B(ω);

  0 ( ) 1/ ( ) MZ     ).
The WAL rotational invariants comprise seven independent (I 1 , I 2 , I 3 , I 4 , I 5 , I 6 and I 7 ) parameters and one dependent (Q) parameter. They can be represented by Mohr circle diagrams [START_REF] Lilley | Magnetotelluric analysis using Mohr circles[END_REF] (Figure 1), and, except I 1 and I 2 , they are taken as sines of angles, which implies an ambiguity in the quadrant to which each angle belongs. Also with the exception of I 1 and I 2 , they are dimensionless and normalized to unity, with their vanishing having a physical interpretation that is related to the geoelectric dimensionality (see [START_REF] Weaver | Characterisation of the magnetotelluric tensor in terms of its invariants[END_REF], for a full description of the invariants).

WAL dimensionality criteria, based on the vanishing or not of some of the invariants (I 3 to I 7 ), are summarized in Table 1. Dimensionality analysis using WAL criteria has been implemented, including data errors and band averages [START_REF] Martí | Geoelectric dimensionality in complex geologic areas: application to the Spanish Betic Chain[END_REF], in the WALDIM code [START_REF] Martí | WALDIM: A code for the dimensionality analysis of magnetotelluric data using the rotational invariants of the magnetotelluric tensor[END_REF]. Given that on field, therefore noisy data, the invariants are rarely precisely zero, the program uses two threshold values (as suggested by [START_REF] Weaver | Characterisation of the magnetotelluric tensor in terms of its invariants[END_REF]:  for I 3 to I 7 ; and  Q , for invariant Q; below which the invariants are taken to be zero.

A c c e p t e d M a n u s c r i p t 6 It is also important to note the parameters that can be derived from the invariants for specific dimensionality cases: In 1D cases, invariants I 1 and I 2 provide information about the 1D magnitude and phase of the geoelectric resistivity ( 1D  and 1D  ). In 2D, the strike angle (referred to as  2D ) can be obtained from the real and imaginary parts of the MT tensor, with  1 and  2 giving the same value for the strike angle (see expressions in the Appendix). In 2D cases affected by galvanic distortion (identified as 3D/2D), the strike angle ( 3D/2D ) is computed considering both the real and imaginary parts of the MT tensor and the distortion parameters, as  1 and  2 [START_REF] Smith | Understanding telluric distortion matrices[END_REF], which are linear combinations of the [START_REF] Groom | Decomposition of the magnetotelluric impedance tensor in the presence of local three-dimensional galvanic distortion[END_REF] twist and shear angles ( It must be remembered that the WAL criteria, as well as the other dimensionality analysis methods, are based on the assumption that the geoelectrical structures are isotropic.

Another tool used to infer the dimensionality in isotropic media is the phase tensor [START_REF] Caldwell | The Magnetotelluric Phase Tensor[END_REF], which is not affected by galvanic distortion (hence only 1D, 2D

and 3D cases can be identified We will compare the phase tensor analysis with the WAL dimensionality criteria for some of the examples presented below.

Electrical anisotropy and modelling

The properties of an anisotropic medium need to be expressed in a tensor form. For the case of electrical anisotropy, the conductivity (σ, reciprocal of the resistivity ρ, σ = 1/ρ) adopts the general form of a symmetric tensor with non-negative components,

xx xy xz yx yy yz zx zy zz                 , (1)
where x (North), y (East) and z (vertically downwards) are the orthogonal axes of a Cartesian coordinate system. The conductivity tensor can represent an intrinsic property of the material (microscopic anisotropy) [START_REF] Negi | Anisotropy in Geoelectromagnetism[END_REF], or it can represent the result of mixing in a preferred orientation of two or more media with differing conductivities (macroanisotropy) (e.g. [START_REF] Wannamaker | Anisotropy versus heterogeneity in continental solid Earth electromagnetic studies: Fundamental response characteristics and implications for physicochemical state[END_REF]. The resolving power of the MT method and the depths at which anisotropic media are typically located (lower crust, upper mantle), usually make it impossible to distinguish between them [START_REF] Weidelt | 3D conductivity models: implications of electrical anisotropy[END_REF].

Using Euler's elementary rotations the conductivity tensor can be diagonalised and its principal directions obtained, namely the strike, dip and slant anisotropy angles (

   ).
In the latter, anisotropy can only be identified by the vertical component of the electric and magnetic fields [START_REF] Negi | Anisotropy in Geoelectromagnetism[END_REF].

For anisotropic media, the MT forward problem must be solved, in general, using a numerical approach. The code of [START_REF] Pek | Finite-difference modelling of magnetotelluric fields in twodimensional anisotropic media[END_REF] uses the finite-difference method to obtain the responses for 1D and 2D anisotropic media.

The magnetotelluric responses obtained from an anisotropic medium are characterized by resistivity shifts, phase splits (which are related to anisotropy contrasts rather than bulk anisotropy of the medium [START_REF] Heise | Anisotropy and phase splits in magnetotellurics[END_REF], and induction arrows not correlated to the principal direction indicated by the MT tensor [START_REF] Pek | Finite-difference modelling of magnetotelluric fields in twodimensional anisotropic media[END_REF][START_REF] Weidelt | 3D conductivity models: implications of electrical anisotropy[END_REF].

Dimensionality analysis of synthetic anisotropic model responses using

WALDIM

In this section we present some examples for synthetic models with anisotropy, the responses of which have been calculated using the code of [START_REF] Pek | Finite-difference modelling of magnetotelluric fields in twodimensional anisotropic media[END_REF] The models were chosen to increase gradually in complexity starting from the most simple. Only 2D situations, not 3D situations, are considered in this study as it is not possible to separate the imprint of anisotropy from 3D effects.

Except when indicated, all the models have dimensions of 860 km (y, towards East) by 186 km (z, vertical downwards), and are discretised using 40 (y) by 30 (z) cells, plus 11 air layers. The responses were computed at each surface node, at the periods indicated in the sections below, following the e +iwt convention for the time-harmonic factor of the electric and magnetic fields. WALDIM analysis was performed for each resulting MT tensor, with 1% random noise having been added to each component. Threshold values of  = 0.1 and  Q = 0.1, which were tested to be consistent with the noise level applied, were used. We also tested, by representing the results using Mohr circle diagrams (following Lilley, 1998b), that the invariant values are obtained as sines of positive angles within the range 0º -90º; and that the dimensionality description obtained from the invariant parameters and Mohr circles are consistent.

ANISOTROPIC HALF-SPACE

For the simplest cases, we considered three models consisting of anisotropic halfspaces. angles has been considered (see Table 2). The responses for each model were computed at T = 1 s, 3.2 s, 10 s, 32 s, 100 s and 320 s.

For the three models, the responses are site independent, and only show slight variations with period due to numerical inaccuracies. Apparent resistivity values depend on the projection of the anisotropy direction on to the x and y axes, as shown in Figure 3a.

Phase values (not shown in the figure) of the off-diagonal components are 45º (xy polarization) and -135º (yx polarization) (as expected from a medium without vertical variations in resistivity). For model 1a, xx and yy apparent resistivities are zero, and hence, the corresponding phases are undetermined. In contrast, for models 1b and 1c, xx and yy phases are 45º and -135º respectively. These responses, observed at a particular site, could be interpreted as a case of galvanic distortion (with shear and anisotropy effects) over a homogeneous medium.

Regarding the dimensionality analysis results, invariant values present similar relationships for the three models, for all sites and frequencies: I 3 = I 4 > , I 5 and I 6 <  and Q <  Q . The exception is I 7 , with values either below or above the threshold value, due to random noise effects. The WAL criteria define the dimensionality as 2D for all models (Figure 3a) and the strike directions are well defined:  1 ≈  2 (=  2D ), with small errors, also due to the noise added. For models 1a and 1b, these angles are coincident The dimensionality and the strike direction also agree with the Mohr diagrams ( hr ≈  hq ≈  1 ≈  2 = 12º), as shown for model 1c in Figure 3b.

For models 1b and 1c, for which the anisotropy directions are not aligned with the measuring axes, two particular features are observed:  3D/2D values (which cannot be represented using Mohr circles) are unstable and are different from  2D (Figure 3a). This does not happen in isotropic 2D structures. In the appendix, the analytical expressions used to obtain the strike directions for the magnetotelluric tensors corresponding to a 2D isotropic model and an anisotropic half-space are developed. In the anisotropic case, the value of  3D/2D is indeterminate, but in the responses of the synthetic model its values are unstable due to the effects of the noise. The main result is that both the analytic expressions and the responses prove that  2D and  3D/2D are not coincident in the case of an anisotropic half-space.

For the three models, phase tensors [START_REF] Caldwell | The Magnetotelluric Phase Tensor[END_REF][START_REF] Heise | Anisotropy and phase splits in magnetotellurics[END_REF] would be represented by unit circles independent of the orientation of the principal directions, and would thus provide no hint of anisotropy.

In model 1a, the fact that all site responses are the same whilst the dimensionality is 2D

indicate that either all the measuring sites are aligned with the strike direction or that the structure is not isotropic but anisotropic. Hence, when the anisotropic directions are coincident with the measuring axes, the responses do not allow the presence of anisotropy in a half-space to be distinguished. In contrast, when anisotropy is not A c c e p t e d M a n u s c r i p t aligned with the measuring axes (models 1b and 1c), the non agreement between the values of the strike directions  2D and  3D/2D is an indication that the half-space over which the measurements are obtained is indistinctly anisotropic. This is an important result, given that it is common to state that 1D anisotropic media are indistinguishable from 2D isotropic media. This type of anisotropic structure cannot be identified using the phase tensor.

1D MEDIA WITH ONE AND TWO ANISOTROPIC LAYERS

The   Ω•m) and α S (45º) were changed in the new layer. The dimensionality pattern for both (Figure 5) is, from the shortest to the longest period: 1D (corresponding to the first isotropic layer), 2D with a 30º strike direction (corresponding to the first anisotropic layer), 3D (due to an abrupt increase in the value of invariant I 7 caused by the inclusion of the second anisotropic layer), and finally 2D, with an approximately 39º strike, a value in between the two anisotropy strike values of the two layers (30º and 45º). In all the 2D cases, as had happened in the case with a single anisotropic layer, the directions  2D and  3D/2D are coincident.

We can summarize that in a 1D medium with one anisotropic layer, dimensionality is 2D with a well defined angle  2D (equivalent to α S , or a projection of the anisotropic directions onto the horizontal if other rotations have been performed), which has the same value as  3D/2D , as would happen in an isotropic medium. In this case the only hint that anisotropy is present is the fact that the responses are the same at all sites, except when the anisotropy angle is 0º, for which responses are equivalent to those of a 2D

A c c e p t e d M a n u s c r i p t 14 model with measurements along the strike direction. When two different anisotropic layers are considered, the dimensionality varies with period: from 2D (corresponding to the first anisotropic layer), to 3D, and back to 2D.

2D ANISOTROPIC MEDIA

In this section, we considered two groups of models based on the examples used in [START_REF] Reddy | Magnetotelluric response of laterally inhomogeneous and anistotropic media[END_REF] and [START_REF] Heise | Anisotropy and phase splits in magnetotellurics[END_REF]. The first group contains models in which the electrical properties vary only in the horizontal direction; the models in the second group possess more general two-dimensional variations.

-Anisotropic dyke:

The models in the first group (Figure 6

( '3 xx   Ω•m, ' 10 yy   Ω•m, ' 20 zz   Ω•m, and º 0 , º 0 , º 30    L D S    )
sandwiched by an isotropic medium of 100   Ω•m. In the third model, both the dyke and the surroundings are anisotropic. The responses for the three models were computed at one period per decade between T = 1 s and T = 10000 s.

For the isotropic model, 3a, the dimensionality is 1D at sites located outside and far from the dyke (Figure 6a). Inside and surrounding the dyke, the dimensionality is 2D (0º strike), except at the first periods for the sites located at the centre of the dyke for which the dimensionality is 1D. At these periods these sites are too far from, and hence not affected by, the dyke boundaries.

A c c e p t e d M a n u s c r i p t

For model 3b (anisotropic dyke surrounded by an isotropic medium), the dimensionality pattern outside the dyke is similar to that of model 3a (Figure 6b): mainly 1D and 2D

(with 0º strike). At the edges of the dyke and at the longest periods the dimensionality is 3D/1D2D. For sites located over the dyke is the dimensionality is 3D, except for the shortest periods at the central part, for which with the dimensionality is 2D with a strike of 30º. In these 2D cases, the strike direction is coincident with the anisotropy angle α S .

However, the direction given by  3D/2D has a different value (60º) from that of  2D (30º), in contrast to what was observed for the anisotropic half-space models (models 1b and 1c), and the distortion parameters are not negligible (φ t = 2º and φ e = -14º).

When both the dyke and surroundings are anisotropic (model 3c), the dimensionality is more complex. Nevertheless there are clear differences with to the results observed for the previous models 3a and 3b, and there are distinctive features associated with each region of the model (Figure 6c). Outside the dyke and far from its edges the dimensionality is 2D with  2D = 55º, which is different from the value given by  3D/2D

(with variable values, as shown in Figure 6c). Still outside but closer to the dyke edges, the dimensionality is mainly 3D/2D with a strike direction of around 75º or 80º and distortion parameters φ t negligible and φ e = -10º. At the edges, is the dimensionality is either 3D/2D or 3D/1D2D. The 3D/2D cases obtained both outside the dyke and at the edges have the peculiarity that, according to the invariant parameters, the dimensionality should be 2D, but the strike directions are inconsistent ( 1 ≠  2 ). It would therefore not be possible to rotate and obtain a regional 2D tensor. Instead, the impedance tensor is better described as 3D/2D with the  3D/2D strike and distortion angles that are small but not negligible [START_REF] Martí | WALDIM: A code for the dimensionality analysis of magnetotelluric data using the rotational invariants of the magnetotelluric tensor[END_REF]. The use of these strike and From the above dimensionality description, the presence of anisotropy can be recognized in the 2D cases, for which the strike directions given by  2D (which agree with the anisotropic azimuth) and  3D/2D are different (Figure 6). Moreover, there are also cases that should be 2D according to the invariants, but for which

 1 ≠  2 .
Therefore, these cases are described as 3D/2D, with the strike direction ( 3D/2D , computed from the real and imaginary parts of the tensor) close to the sum of both anisotropic directions (80º).  hr and  hq do not provide the correct strike direction either, as they are computed using the real or the imaginary parts separately. Also, for sites over the dyke and at the edges, some 3D/1D2D and 3D cases are obtained.

-2D conductive bodies and an anisotropic layer:

The second set of models is taken from the 2D examples used in [START_REF] Heise | Anisotropy and phase splits in magnetotellurics[END_REF].

This set explores phase splits in responses from anisotropic structures, and the identification of anisotropy using phase information (in particular, the phase tensor).

Model 4a contains an anisotropic layer (with main directions along the x, y and z axes)

and two conductive blocks. .Model 4b is isotropic and contains two conductive blocks similar to those in model 4a (Figure 7). [START_REF] Heise | Anisotropy and phase splits in magnetotellurics[END_REF] show how both models give Rotational invariants and dimensionality of the responses of these two models were computed between 1 s and 30000 s (Figure 8). The invariants have similar values for both models, and hence the dimensionality displays similar patterns for both models. In general, the dimensionality is 1D for periods up to approximately 100 s, and 2D for the rest. However, the dimensionality of the first and last sites, located on top of the conductive blocks, is different for each model. For model 4a dimensionality is 1D up to 100 s, then becoming 2D as a consequence of the directionality introduced by the anisotropic layer, affecting all sites. For model 4b all the cases are 1D as these sites are not affected by the lateral contrasts at the limits of the two conductive bodies. The phase tensor ellipses, which were computed for both models for the first site, also show this difference between the two models (Figure 9). These results confirm that for 2D models with anisotropic structures aligned with the main directions, both the invariants and the phase tensor provide the same information, and that for this example they cannot distinguish between the anisotropic and isotropic models.

Additionally, we considered model 4a and modified the anisotropic layer by applying a rotation of the principal directions (α S = 30º). The dimensionality of the responses of the resulting model, identified as 4c, is shown in Figure 8c. The dimensionality pattern is significantly more complex than for the previous models. Up to 100 s, the dimensionality is similar to that of models 4a and 4b (mostly 1D with 2D cases at the rightmost side of the model due to the shallow conductive structure). For periods around (with an approximately 15º strike). At longer periods, the general trend is that the cases that were 1D and 2D in models 4a and 4b become 2D and 3D, respectively; with some 3D/2D and 3D/1D2D exceptions. In all 3D/2D cases (most of them at 100 s), as happened for the model with the anisotropic dyke (3c), invariant values indicate 2D

dimensionality, but, given that the two strike directions,  1 and  2 , are significantly different, the impedance tensor is better described as 3D/2D with  3D/2D . This observation is a clear indication of the presence of anisotropy in the structures, with anisotropic directions non-aligned to the principal structural directions. In the phase tensor diagrams of these model responses (Figure 10) an equivalent effect can be observed at 100 seconds: the values of  are negligible (note that only angle values lower than 3º are considered negligible), whereas the main directions of the ellipse differ significantly from the strike angle .

Hence, for 2D models, both the WALDIM criteria and the phase tensor diagrams are able to identify the presence of anisotropic structures with principal directions not coincident with the measuring axes.

Anisotropy in field data: the COPROD dataset

In this final section we refer to one case of field data that has been associated with anisotropy. This is the well known COPROD2 dataset, from southern Saskatchewan and Manitoba (Canada), which revealed the presence of the North American Central Plains conductivity anomaly (NACP) (Jones andCraven, 1990, Jones et al., 1993). This dataset was used to test inversion codes (see [START_REF] Jones | The COPROD2 dataset: Tectonic setting, recorded MT data and comparison of models[END_REF]. Some of the 2D models that We computed the dimensionality for the synthetic tensors of the sites located over the anisotropic body. We obtained 1D cases, and 2D cases with 0º strike (anisotropy aligned with the measuring directions) (Figure 11a).

The data for site 85_314 was used by [START_REF] Martí | WALDIM: A code for the dimensionality analysis of magnetotelluric data using the rotational invariants of the magnetotelluric tensor[END_REF] as an example for dimensionality analysis using the WALDIM code. Up to 10 s, the data can be described as 1D. At periods longer than 10 s invariant values indicate 2D. However, at these periods strike angles  1 and  2 differ significantly, and hence the data were better described as 3D/2D, with a strike direction around 80º and small twist and shear distortion angles (lower than 5º) (Figure 11b). This allowed 2D regional tensors to be obtained from tensor decomposition. According to the tests presented here, the discrepancy between the dimensionality descriptions from the model with the anisotropic block and the field data lies in the fact that in the synthetic data all the diagonal responses are null, whereas in the field data the values of the diagonal components, especially for the longest periods, are significant. From our new characterization of anisotropy in dimensionality analysis, it is clear that the dimensionality of site 85_314 is compatible with a 2D model that contains at least one anisotropic block or layer, having anisotropy directions aligned with the strike indicated in the dimensionality analysis (in this case of around 80º). Hence, if the anisotropic block modeled by [START_REF] Jones | Electromagnetic interrogation of the anisotropic Earth: Looking into the Earth with polarized spectacles[END_REF] had an anisotropic azimuth of 80º, the invariant values of the responses would correspond to 2D structures with two different strike directions  1 and  2 . This would be in agreement with the observed data.

WAL criteria extended to accommodate anisotropy

The results obtained from this study have allowed specific relationships to be established between the invariants and strike directions that are linked to the presence of anisotropy. In general, these conditions are not recognized from a single tensor alone, but from the pattern at different sites and periods. The main imprint of anisotropy can be seen in the 2D cases (according to WAL isotropic criteria), with strike directions that are not consistent, or relationships between tensors that would not correspond to isotropic structures. In these cases, the strike obtained is related to the orientation of the anisotropy rather than to the structural direction. Table 3 contains the new dimensionality criteria extended to accommodate these cases with anisotropy and to distinguish them from isotropic two-dimensionality: anisotropic half-space, a 1D medium with one anisotropic layer, and an anisotropic 2D medium.

However, it must be remembered that it is not always possible to identify anisotropy when the main directions are aligned with the measuring axes, or to retrieve all the Table 3 considers the dimensionality observed in a particular tensor. In particular situations described in the text, some patterns can be observed such as that of a 1D model with two anisotropic layers (2D, 3D and 2D cases, as the period increases).

Hence, once the dimensionality of the full dataset is obtained (it is recommended to plot dimensionality maps), one should check for anisotropic imprints and patterns as described in the text, and evaluate what type of anisotropic media might exist beneath the survey area.

Conclusions

The most important contribution of this study is the demonstration that it is possible to identify the presence of anisotropy in the dimensionality description given by the WAL criteria. In addition, we have extended the WAL invariants criteria to differentiate anisotropic from isotropic media. Hence, when assessing the dimensionality of a dataset that is considered to contain anisotropy effects, one should follow the original WAL criteria (Table 1), plus the new conditions described in Table 3. The exception is when the principal anisotropy directions are aligned with the measuring axes. In this situation, if the anisotropic media is 2D, the information contained in the induction arrows might be useful.

A c c e p t e d M a n u s c r i p t

22 Another important point is that, except in very simple cases, the anisotropy cannot be identified from one site alone. It is fundamental to check for the consistency of dimensionality with neighbouring sites or periods.

Finally, the comparison of the dimensionality description obtained using the WAL invariant criteria with that from phase tensor diagrams allowed us to conclude that, in some cases, both provide the same information. However, when the phases do not change with period, such as in the case of an anisotropic half-space, only the WAL criteria enable the anisotropy to be identified. It is also important to note that in some cases the strike angle can only be computed from  3D/2D , which considers the real and imaginary parts of the tensor, as opposed to the direction defined from the Mohr circles,  h , which uses the real or the imaginary parts separately.

Acknowledgements

The authors thank Ted Lilley for his critical review, which greatly helped improving the manuscript. We acknowledge Josef Pek for providing the 2D anisotropic forward code [START_REF] Pek | Finite-difference modelling of magnetotelluric fields in twodimensional anisotropic media[END_REF]. This work has been funded by projects CGL2006-10166 and CGL2009-07604. A. Martí thanks the Universitat de Barcelona and the Department of Earth Sciences at Memorial University of Newfoundland (MUN) for facilitating her research term at MUN.

APPENDIX A

In this appendix we first summarize the expressions used to compute the strike directions from the magnetotelluric tensor using [START_REF] Weaver | Characterisation of the magnetotelluric tensor in terms of its invariants[END_REF] notation.

A c c e p t e d M a n u s c r i p t

23 Secondly, we derive these expressions for the theoretical magnetotelluric tensors corresponding to A) a 2D isotropic structure, rotated an angle θ from the strike direction, and B) an anisotropic half-space, with the main anisotropic directions rotated an angle α S .

1. Strike expressions:

The complex parameters

MM M i MM                                                                 . (A1)
If the tensor corresponds to a 2D structure, the strike direction (θ 2D ) can be computed from using either the real or the imaginary parts of 2  and 3  , which lead to the same result:

θ 2D = θ 1 = θ 2 :   3 1 2 tan 2     , (A2) and   3 2 2 tan 2     . (A3)
A c c e p t e d M a n u s c r i p t

Both the i  parameters and the angles θ 1 and θ 2 can be represented in Mohr circle diagrams (for the real and imaginary parts), which are also used to represent WAL invariants (Figure 1).

If the 2D structure is affected by galvanic distortion, the strike direction (θ 3D/2D ) can be recovered using the expression:

  12 34 3 / 2 13 24 tan 2 DD dd dd     , ( A4
)
where 12 i j j i ij d II      
, and I 1 and I 2 are rotational invariants of the MT tensor.

Given that 2D is a particular case of 3D/2D (where the galvanic matrix is the identity), the same expression works to compute the strike, so that:

θ 2D = θ 1 = θ 2 = θ 3D/2D .

Particular cases:

A. 2D isotropic structure:

2 0 0 xy D yx M M M        , ( A5
)
if the tensor is rotated an angle θ:

22 2 22 ( )sin •cos •cos •sin ' • • •sin •cos ( )•sin •cos xy yx xy yx T D xy yx xy yx M M M M M R M R M M M M                    , ( A6
M M M M M M M M d d M M M M M M M M                   2 2 2•sin cos sin(2 ) tan(2 ) sin cos cos(2 )              
, which proves that: The analytic expression of the MT tensor corresponding to an anisotropic half-space, with the main anisotropic directions rotated an angle α S is obtained using the development from [START_REF] Pek | Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic media[END_REF]:

3 / 2 1 2 D D        .
  2 2 1 2 2 S S anis S S d•sin( ) s d•cos( ) M C• • i s d•cos( ) d•sin( )                ,
where C is a constant,

xx yy s ' '     and xx yy d ' '     . 1 2 3 4 0 • cos(2 )(1 ) 2 • sin(2 )(1 ) 2 • •(1 ) S S d C i d C i C s i              
Given that both the real and imaginary parts have the same value,

1 2 sin(2 ) 2 tan(2 ) tan(2 ) tan(2 ) cos(2 ) 2 S S S d d           , which proves that 2 1 2 D S        .
On the other hand, if the strike direction is computed using the expression of θ 3D/2D , d ij = 0, for any i, j because real and imaginary parts of the tensor are identical. Consequently: by [START_REF] Jones | Electromagnetic interrogation of the anisotropic Earth: Looking into the Earth with polarized spectacles[END_REF], which fits the off-diagonal components of site 85_314 in the COPROD2 dataset; and b) for all the components of site 85_314 from the COPROD2 dataset (modified from [START_REF] Martí | WALDIM: A code for the dimensionality analysis of magnetotelluric data using the rotational invariants of the magnetotelluric tensor[END_REF].

Table captions:

Table 1: Dimensionality criteria according to the WAL invariant values of the magnetotelluric tensor (modified from [START_REF] Weaver | Characterisation of the magnetotelluric tensor in terms of its invariants[END_REF]. Row "2D" with the grey background is extended in Table 3, where structures with anisotropy are considered.

Table 2: Resistivity values and orientations of the three anisotropic half-space models 1a, 1b and 1c.

  de Geodinàmica i Geofísica. Universitat de Barcelona, Spain. 2 : Earth Sciences Department. Memorial University of Newfoundland, Canada. annamarti@ub.edu Submitted to Physics of the Earth and Planetary Interiors Version 2. July 2010

  presented a Mohr circles analysis for data with phases out of quadrant, although not particularly related with anisotropy.

  diagrams to represent the responses of models with electric anisotropy.



  around x'-axis and L  around z''-axis) (Figure2). Hence, the anisotropy, specified in terms of the relationships between the components along the principal directions of the conductivity tensor, are azimuthal anisotropy (''

  In

  have performed the dimensionality analysis using the WALDIM code and we have analyzed the results indicating which features are characteristic of the anisotropic structures.

  , and are distinguished from each other by the orientation of the principal directions. In the first model (1a), these coincide with the measurement axes. In the second (1b), these have been rotated through a strike angle 40 S   around the z axis. In the third model (1c), a general rotation using dip (55º) and slant (20º)

  (0º deg or 40º respectively), and for model 1c it is 12º, due to the projection onto the horizontal plane of the new x' and y' axes, resulting from the dip and slant rotations.

  first one-dimensional model presented here (model 2a) was taken from one of the examples provided with the Pek and Verner (1997) code. It consists of a layered structure with an embedded anisotropic layer (Figure 4): ( '1 xx  model responses were computed at 10 periods between T = 1 s to T = 32000 s. The MT responses, which are shown in Figure 4, are the same at all sites.. Diagonal responses are coincident (xx = yy), whereas the off-diagonal responses show a split between the polarizations. The off-diagonal resistivity and phases are plotted together with the responses (xy = yx) of two models in which the anisotropic layer is replaced with an isotropic one; the first model with a 1 Ω•m layer, and the second with a 100 Ω•m layer (Figure 4). Because of the rotation (α S ) of the principal directions, the values of the off-diagonal resistivities and phases for the model with the anisotropic layer are smoother than those for the models with the isotropic layers. The WAL dimensionality criteria (I 3 = I 4 >  , I 5 , I 6 and I 7 <  and Q >  Q ) indicate 2D structures with  2D = 30º (= α S ) for all periods (Figure 4), except for T=1 s at which the A c c e p t e d M a n u s c r i p t 13 criteria indicate 1D structure because the skin depth (5 km) is smaller than the top of the anisotropic layer. For the periods at which 2D structure is indicated, the strike direction computed as  3D/2D is coincident with  2D and the distortion parameters are practically null. The effects of the inclusion of a second anisotropic layer just below the first one were also investigated by considering the third layer of model 2a to be anisotropic as well. In the first of these models (model 2b), this new anisotropic layer has the same resistivity values as the upper one, but with the main directions rotated at an angle α S = 45º. In the second model (model 2c), both the resistivity values ( '1 xx   Ω•m and ' ' 10 yy zz

  ) consist of a vertical dyke intruded into a medium with differing electrical properties. Initially we consider a model in which both the dyke and the surroundings are isotropic (model 3a, ). A second model (model 3b) consists of an anisotropic dyke

  , in isotropic structures, the decomposition of the impedance tensor to be performed and a 2D tensor recovered. Inside the dyke, at the shortest period, the dimensionality is 2D, with  2D = 30º = α S , inconsistent with  3D/2D (70º), and with non-negligible values of the shear distortion angle (φ e = -10º). As the period increases, the dimensionality becomes 3D, 3D/2D (with  3D/2D = 75º, φ t = -12º and φ e = -10º), and finally 3D/1D2D.

  and induction arrow responses, except at the longest periods, where the induction arrows for the isotropic model are significant, whereas in the anisotropic model they are almost null.

  most of the 1D cases become 2D ( 2D being consistent with  3D/2D ) or 3D/2D

  of multiple isotropic high conductivity bodies separated by resistive regions. Jones (2006) revisited the data and, using one of the sites on top of the NACP anomaly as a reference (85_314), proposed a 2D anisotropic model. This model consists of a thin superficial conductive layer (3 Ω•m), a 100 km thick lithosphere of 1000 Ω•m, in which a single anisotropic block ( ) is embedded, and a basal conducting layer of 10 Ω•m. The offdiagonal responses for this model are in good agreement with those of the observed data, reproducing the split between TE and TM modes.

  anisotropy from the observed responses and the dimensionality analysis alone.

Figure 1 :

 1 Figure captions:

Figure 2 :

 2 Figure 2: Diagram of successive Euler rotations applied to generate any orientation of

Figure 3

 3 Figure 3: a) Dimensionality and apparent resistivity responses for the three anisotropic

Figure 4 :

 4 Figure 4: Cross section of model 2a, corresponding to a layered model with an

Figure 5 :

 5 Figure 5: Top: Cross sections of models 2b and 2c, consisting of 1D models with two

Figure 6 :

 6 Figure 6: Cross-section of models 3a, 3b and 3c and the corresponding dimensionality

Figure 7 :

 7 Figure 7: Cross-sections of models 4a and 4b (from Heise et al., 2006), used to compute

Figure 8 :

 8 Figure 8: Dimensionality patterns corresponding to the responses of models 4a and 4b

Figure 9 :

 9 Figure 9: Phase tensor diagrams corresponding to site 1 (located at km 0) for models 4a

Figure 10 :

 10 Figure 10: Phase tensor diagrams corresponding to the responses of model 4c, for 5

Figure 11 :

 11 Figure 11: Dimensionality cases for: a) the responses of the anisotropic model presented

  Figure 1 Click here to download high resolution image

  Figure 2 Click here to download high resolution image

  Figure 3 Click here to download high resolution image

  Figure 4 Click here to download high resolution image

  Figure 5 Click here to download high resolution image

  Figure 9 Click here to download high resolution image

(

  Figure 10 Click here to download high resolution image

  •

								j    j	i		j	(j = 1,4), are defined as linear combinations of
	the magnetotelluric tensor components:	   1 xx MM	yy		/2	,	   xy 2 MM	yx		/2	,
	   3 xx MM	yy		/2	and	4  		xy MM 	yx		/2	:
	xx	xy			1	3	2		4			1		3	2	4	1	3	2	4
	yx	yy			2	4	1		3			2		4	1	3	2	4	1	3

Table 3 :

 3 Dimensionality criteria extended to anisotropic structures, characterized by the WAL invariants criteria indicating isotropic 2D.

Table 1 .

 1 Martí et al.

Table 1 A

 1

	'   xx	Homogeneous models 50 Ω•m and ' ' 500 yy zz  	Ω•m		Anisotropy angles
		Model 1a		S		0º ,		D		0º ,		L		0º
		Model 1b		S			40º ,		D		0º ,		L		0º
		Model 1c		S		0º ,		D			55º ,		L		20º
			M a n u s c r i p t
		c c e p t e d

Table 2 .

 2 Martí et al.

Table 2

 2

Table 3 . Martí et al.

 3

Table 3

 3