F Paglialunga

F Passelègue

M Lebihain

M Violay

Frictional weakening leads to unconventional singularities during dynamic rupture propagation

Earthquakes i.e. frictional ruptures, are commonly described by singular solutions of shear crack motions. These solutions assume a square root singularity order around the rupture tip and a constant shear stress value behind it, implying scale-independent edge-localized energy. However, recent observations of large-scale thermal weakening accompanied by decreasing shear stress potentially affecting the singularity order can challenge this assumption. In this study, we replicate earthquakes in a laboratory setting by conducting stick-slip experiments on PMMA samples under normal stress ranging from 1 to 4 MPa. Strain gauges rosettes, located near the frictional interface, are used to analyze each rupture event, enabling the investigation of shear stress evolution, slip velocity, and material displacement as a function of distance from the rupture tip. Our analysis of the rupture dynamics provides compelling experimental evidence of frictional rupture driven by enhanced thermal weakening. The observed rupture fronts exhibit unconventional singularity orders and display slip-dependent breakdown work (on-fault dissipated energy). Moreover, these findings elucidate the challenges associated

Introduction

Frictional rupture phenomena, including natural earthquakes, are often described by singular solutions of shear crack motions [START_REF] Freund | The mechanics of dynamic shear crack propagation[END_REF][START_REF] Palmer | The growth of slip surfaces in the progressive failure of over-consolidated clay[END_REF][START_REF] Rice | The mechanics of earthquake rupture[END_REF]. For such cracks, the stress field at the rupture tip is described by a square root singularity (ξ = -0.5), constant residual stress is expected far behind the rupture tip, and the energy balance condition equates the energy release rate G (i.e. rupture growth driving force) to a constant value of fracture energy G c (i.e. resistance to rupture growth).

This was confirmed by experimental and numerical observations, where the onset of frictional sliding, the evolution of the rupture speed, and the rupture length were predicted by Linear Elastic Fracture Mechanics (LEFM) [START_REF] Bayart | Fracture mechanics determine the lengths of interface ruptures that mediate frictional motion[END_REF][START_REF] Kammer | Linear elastic fracture mechanics predicts the propagation distance of frictional slip[END_REF][START_REF] Kammer | Fracture energy estimates from largescale laboratory earthquakes[END_REF][START_REF] Svetlizky | Classical shear cracks drive the onset of dry frictional motion[END_REF][START_REF] Xu | Robust Estimation of Rupture 29 Properties at Propagating Front of Laboratory Earthquakes[END_REF], and suggesting that the fracture energy controlling the dynamics of the rupture tip might be an interface property.

Such an analysis often relies on the hypothesis of negligible frictional weakening far away from the rupture tip (i.e. outside of the cohesive zone). However, it is widely recognized that fault shear stress is likely to evolve during seismic slip due to (i) velocity and slip dependencies [START_REF] Marone | Laboratory-derived friction laws and their application to seismic faulting[END_REF], (ii) activation of thermal weakening processes [START_REF] Toro | Fault lubrication during earthquakes[END_REF][START_REF] Hirose | Growth of molten zone as a mechanism of slip weakening of simulated faults in gabbro during frictional melting[END_REF][START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF], (iii) dilatancy inducing fluid pressure changes [START_REF] Brantut | Dilatancy-induced fluid pressure drop during dynamic rupture: Direct experimental evidence and consequences for earthquake dynamics[END_REF][START_REF] Rice | Earthquake precursory effects due to pore 28 fluid stabilization of a weakening fault zone[END_REF][START_REF] Segall | Dilatant strengthening as a mechanism for slow slip events[END_REF].

These changes in the residual stress behind the rupture tip could induce a slip dependency of the apparent fracture energy (nowadays more commonly called breakdown work [START_REF] Tinti | Earthquake fracture energy inferred from kinematic rupture models on extended faults[END_REF]) estimated for natural earthquakes [START_REF] Abercrombie | Can observations of earthquake scaling constrain slip weakening?[END_REF][START_REF] Lambert | Rupture-dependent breakdown energy in fault models with thermo-hydro-mechanical processes[END_REF], in contrast to the LEFM definition. The breakdown work W bd is a quantity commonly used to study the energy balance of earthquakes and is defined as an energy term including all on-fault dissipative processes W bd = D 0 τ (D ′ ) -τ min dD ′ , with τ the shear stress acting on the fault, τ min the minimum shear stress reached on-fault, and D the fault slip. It can be observed that, by definition, W bd is a slip-dependent quantity. It is therefore important to be aware of how possible stress weakening may affect rupture dynamics and the energy release that controls it. Such weakening phenomenon has been observed in high-velocity friction experiments (for a comprehensive review, see [START_REF] Toro | Fault lubrication during earthquakes[END_REF]). Nonetheless, friction experiments seldom exhibit elastodynamic rupture propagation due to the imposed fault slip by the experimental procedure [START_REF] Chen | Dynamic fault weakening during earthquakes: Rupture or friction?[END_REF]. Consequently, the experimental documentation of the coupling between fault weakening and elastodynamic rupture response remains limited [START_REF] Rosakis | Recent Milestones in Unraveling the Full-Field Structure of Dynamic Shear Cracks and Fault Ruptures in Real-Time : From Photoelasticity to Ultrahigh-Speed Digital Image Correlation[END_REF][START_REF] Paglialunga | On the scale dependence in the dynamics of frictional rupture : Constant fracture energy versus size-dependent breakdown work[END_REF].

In these regards, our recent work highlighted that a long-tailed weakening can emerge after a first rapid weakening during frictional rupture experiments [START_REF] Paglialunga | On the scale dependence in the dynamics of frictional rupture : Constant fracture energy versus size-dependent breakdown work[END_REF], resulting in a slip-dependent breakdown work. Despite this observation, the rupture dynamics, analyzed through LEFM, showed to be controlled by a constant fracture energy G c , in agreement with previous studies [START_REF] Bayart | Fracture mechanics determine the lengths of interface ruptures that mediate frictional motion[END_REF][START_REF] Kammer | Linear elastic fracture mechanics predicts the propagation distance of frictional slip[END_REF][START_REF] Kammer | Fracture energy estimates from largescale laboratory earthquakes[END_REF][START_REF] Svetlizky | Classical shear cracks drive the onset of dry frictional motion[END_REF][START_REF] Xu | Robust Estimation of Rupture 29 Properties at Propagating Front of Laboratory Earthquakes[END_REF]. However, analyzing such frictional ruptures in the framework of LEFM relies on the assumption of constant residual stress behind the rupture tip. The observed long-tailed weakening could call into question this assumption and limit the framework's applicability to fully describe frictional ruptures, explaining the observed mismatch between G c and W bd [START_REF] Paglialunga | On the scale dependence in the dynamics of frictional rupture : Constant fracture energy versus size-dependent breakdown work[END_REF].

Moreover, theoretical studies have shown that continuous stress weakening can modify the singularity order controlling the stress and displacement fields around the rupture tip, deviating from the square-root singularity commonly adopted in LEFM, and leading to an unconventional singularity order (ξ ̸ = -0.5) [START_REF] Garagash | Multiscale tip asymptotics in hydraulic fracture with leak-off[END_REF][START_REF] Viesca | Ubiquitous weakening of faults due to thermal pressurization[END_REF][START_REF] Brantut | The fracture energy of ruptures driven by flash heating[END_REF]Brener and Bouchbinder, 2021b). In particular, when frictional ruptures are described by ξ ̸ = -0.5, the stress σ ij and displacement u i fields obey respectively the following scaling relationships (Brener and Bouchbinder, 2021a):

σ ij ∝ K (ξ) r ξ and u i ∝ K (ξ) r (ξ+1) /µ, with K (ξ)
the ξ-generalized dynamic stress intensity factor, r = x -x tip the distance from the rupture tip of a point of observation located on the fault plane, and µ the dynamic shear modulus. These lead to the following relation:

W bd ∝ [K (ξ)
] 2 r (1+2ξ) /µ, valid for r > x c , with x c the cohesive zone size (Eq. ( 5) from (Brener and Bouchbinder, 2021a)). From this relation, it can be easily noticed that for ξ = -0.5, the W bd dependence on r completely vanishes, making the breakdown work independent of the distance from the rupture tip. This does not happen when ξ ̸ = -0.5, for which W bd has a direct dependence on r.

So far, the occurrence of such unconventional singularities during fric-tional ruptures has not been measured at the laboratory scale. In this paper, we present new data analyzed in a recently-derived theoretical framework, demonstrating the first experimental evidence of strain and stress perturbation caused by unconventional singularities associated with velocitydependent frictional weakening. These experimental findings are supported by theoretical explanations about the emergence of unconventional singular fields during dynamic rupture.

Methods

We performed stick-slip experiments in a biaxial apparatus working in a 2D single shear configuration under an applied normal stress ranging from 1 to 4 MPa (Figure 1a). The experimental setup is the same one used and described in [START_REF] Paglialunga | On the scale dependence in the dynamics of frictional rupture : Constant fracture energy versus size-dependent breakdown work[END_REF]. The tested samples consist of two polymethylmethacrylate (PMMA) blocks of dimensions (20x10x3) cm (top block) and (50x10x3) cm (bottom block), generating, once put into contact, an artificial fault of (20x3) cm. The external loading is imposed using two hydraulic pumps. The normal load is applied to the top block and kept constant while the shear load is manually increased and applied to the bottom block inducing, once reached the fault strength, stick-slip events. Strain gages rosettes (oriented along 45 • , 90 • , 135 • to the fault plane), located 1 mm away from the frictional interface, were used to compute the local strain and stress tensors. Note that at such distance from the fault plane, the dynamic propagation of the rupture will induce a stress perturbation (rise of stress ahead of the rupture tip, followed by a stress drop behind it). Such a stress perturbation will be maximum on the fault plane (y = 0+) and decrease away from it (y > 0), as predicted by fracture mechanics. This feature is not observable through macroscopic measurements of strain or stress, hence the need for local ones. The strain tensor rotation was obtained through the conversion of ε 1 , ε 2 , ε 3 into ε xx , ε xy , ε yy following:

ε xy = ε 3 -ε 2 2 , ( 1
)
ε yy = ε 1 , (2)
ε xx = ε 3 + ε 2 -ε 1 (3)
The local strain temporal evolution shows clear perturbations concurrent with stick-slips (Figure 1b). By zooming-in in time, details of the instability can be caught (Figure 1c), showing a first (main) rupture front, followed by a series of secondary fronts probably caused by rupture reflections at the fault edges. To study the rupture dynamics, only the main front was considered, selecting a time window around the first strain perturbation (Figure 1d-e).

Together with the vicinity of the measurement location to the fault, the measuring interval is of paramount importance for properly addressing the elastodynamics of the studied events. From here the need for a sufficiently high recording frequency (2 MHz in this case) and a well-defined and short time interval excluding possible subsequent events. Note that the following analysis and discussions will exclusively focus on the dynamics of the main rupture front for each stick-slip event. All the experimental curves that will be shown will refer to a defined time window, systematically smaller than the expected propagation time along the fault interface (for example, the rupture showed in Figure 1e is described by a temporal window of ∼ 45 µs). The rupture propagation velocity C f was estimated by computing the ratio between the distance among the strain gauge locations and the rupture front travel time from one location to the other. For each event, the particle velocity was then computed through the strain component parallel to the slip direction as ux = -C f ε xx . This estimate has been shown to be comparable to distinct measurements of slip motions associated with the propagation of the seismic rupture in previous experimental studies [START_REF] Svetlizky | Classical shear cracks drive the onset of dry frictional motion[END_REF][START_REF] Paglialunga | On the scale dependence in the dynamics of frictional rupture : Constant fracture energy versus size-dependent breakdown work[END_REF]. The fault slip velocity was considered equal to twice the particle velocity measured through the strain gauges (V = 2 ux ), assuming an anti-symmetrical distribution of slip and slip rate. This assumption seems legitimate given that the two samples have comparable dimensions, the same width, and are made of the same material. Integrating V during the propagation time, local horizontal displacement u x could be estimated as well. The slip displacement D of the fault is computed as twice (refer to the assumption described just above) the horizontal displacement (D = 2u x ), assuming the material displacement measured through the strain gauge 1 mm away from the fault is comparable to the one occurring on-fault. The acquired temporal evolutions of strain were converted into spatial evolutions following the methods used in [START_REF] Svetlizky | Classical shear cracks drive the onset of dry frictional motion[END_REF], and by assuming that the rupture velocity obtained was locally constant during rupture propagation. Moreover, by assuming plane stress conditions, the stress tensor was computed through the elastic properties of PMMA. Given the viscoelastic nature of PMMA, dynamic values of its elastic modulus were used in this study, E = 5.7 GPa, and a Poisson ratio ν = 0.33, both obtained through bench measurements of seismic velocities (C P = 2705 m/s and C S = 1340 m/s).

Results

Each rupture event was studied through the evolution of shear stress, slip velocity, and material displacement as a function of the distance from the rupture tip (Fig. 2). In such spontaneous ruptures, the slip (or slip rate) history results from the coupling of elastodynamics and the constitutive friction law, rather than being prescribed. In all the studied events, local shear stress evolution exhibited an increase ahead of the rupture tip followed by a first significant decrease within the first micrometers of slip and a second mild one within larger distances (Fig. 2a) as recently observed [START_REF] Paglialunga | On the scale dependence in the dynamics of frictional rupture : Constant fracture energy versus size-dependent breakdown work[END_REF]. A rapid increase of slip velocity was observed concurrent with the passage of the rupture front, followed by a slow decay occurring with distance from the rupture tip. The peak slip velocity V max showed a clear dependence with estimated rupture speed, with ∼ 0.08 m/s for C f ≈ 220 m/s up to ∼ 0.8 m/s for C f ≈ 840 m/s (Fig. 2b). The evolution of horizontal displacement u x presented values close to 0µm ahead of the rupture tip (values slightly deviate from 0 due to off-fault measurement) and a sharp increase behind it (Fig. 2c), with final displacements ranging between 3.9 and 28 µm. Subsequently, the fault strength weakening was analyzed through the evolution of the local shear stress τ with the fault's slip displacement D. The fault's weakening presents a sharp decrease of shear stress occurring within the first microns of slip, followed by a milder decrease occurring within a larger amount of slip (Fig. 3a). The breakdown work evolution was computed as where D((x -x tip ) = 0) is the displacement at the passage of the rupture tip.

W bd = D D((x-x tip )=0) (τ -τ (D)) dD (
Since no slip is expected to occur ahead of the rupture tip on the fault plane (C f = 0 when (x -x tip ) > 0), the breakdown work evolution was computed only from slip occurring after the passage of the rupture tip (x -x tip ) = 0, neglecting fictitious contributions due to elastic strain of the bulk at the measurement location. The evolution of W bd showed a first increase with slip described by a slope close to 1 : 2 and a subsequent increase described by a slope of ∼ 1 : 0.6(±0.1) (Fig. 3b), suggesting the existence of anomalous singularities (ξ ̸ = -0.5). Given our interest in studying the second branch of this evolution, the slope of it was measured. This was done by fitting the evolution of W bd with D for D > D c with a first-degree polynomial in log-log scale. Then, ξ was derived from estimates of the power law exponent For each curve, the slope of the second branch of the evolution corresponds to a given value of singularity order. c. Evolution of estimated ξ values with peak slip velocity V max .

through:

W bd (D) = G c D D c ( 1+2ξ 1+ξ ) , (5)
for large slip displacement D ≳ D c . The values were found to range between -0.4 and -0.2 (Fig. 3c). Another possible route proposed by Brener and Bouchbinder (2021a) is to estimate the deviation ∆ξ = ξ + 1/2 of xi from the LEFM value -1/2 from the dependence of the breakdown energy in the distance r = x -x tip from the crack tip. However, this requires prior knowledge of the fracture energy G c and the process zone size x c , which are difficult to constrain from experimental measurements.

4. Theoretical modeling of the kinematic fields around the rupture tip for unconventional singularity order

While the first increase of breakdown work with slip can be explained by a slip-weakening behavior of the fault, the subsequent increase (power law of 1:0.6) is unexpected from the conventional theory of LEFM. If such a continuous weakening stage controlled the dynamics of the rupture, stress fields with a scaling σ ∝ r ξ should be observed behind the rupture tip, as expected from theoretical studies [START_REF] Brantut | The fracture energy of ruptures driven by flash heating[END_REF]Brener and Bouchbinder, 2021b;[START_REF] Garagash | Multiscale tip asymptotics in hydraulic fracture with leak-off[END_REF][START_REF] Viesca | Ubiquitous weakening of faults due to thermal pressurization[END_REF], with the singularity order ξ different from the square root singularity. To further investigate the dynamics of rupture, the temporal evolution of the strain perturbations generated by the passage of the rupture front (∆ε xy , ∆ε xx ) was compared to the theoretical predictions obtained considering both a square root singularity (LEFM) and an unconventional singularity (Brener and Bouchbinder, 2021a).

For the LEFM theoretical prediction, the stress field perturbation around the rupture tip takes the following general form (for a detailed description please refer to [START_REF] Freund | Dynamic fracture mechanics[END_REF]; [START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF]):

∆σ ij (r, θ) = K II √ 2πr Σ II ij (θ, C f ) (6)
where K II the stress intensity factor, and Σ II ij (θ, C f ) the angular variation function. Coordinates are expressed in the polar system with (r, θ) respectively the distance from the crack tip and the angle to the crack's plane.

In the unconventional theory framework, the stress fields were derived from the elastodynamic equations assuming a steady-state rupture velocity.

The equations obtained present the following form:

σ xx (r, θ) = 4(ξ + 1)K (ξ) II √ 2πR(C f ) [α s (1 + 2α 2 d -α 2 s )r ξ d sin(ξθ d ) -α s (1 + α 2 s )r ξ s sin(ξθ s )], (7)
τ (r, θ) = σ xy (r, θ) = τ min + 2(ξ + 1)K (ξ) II √ 2πR(C f ) [4α s α d r ξ d cos(ξθ d )-(1+α 2 s ) 2 r ξ s cos(ξθ s )], (8)
σ yy (r, θ) = - 4(ξ + 1)α s (1 + α 2 s )K (ξ) II √ 2πR(C f ) [r ξ d sin(ξθ d ) -r ξ s sin(ξθ s )], (9)
where and

K (ξ) II = lim r→0 (2 √ 2π) (ξ+1) r -ξ τ (r, 0 +-) is the ξ-generalized stress inten- sity factor, α d = 1-C f C d 2 and α s = 1-C f
r d = r 1 -C f sin(θ) C d 2 , r s = r 1 -C f sin(θ) Cs 2 .
The displacement field related to the unconventional rupture phenomenon can be predicted by (Brener and Bouchbinder, 2021a):

u x (r, θ) = 2K (ξ) II µ √ 2π 2α s r ξ+1 d sin[(ξ + 1)θ d ] -α s (1 + α 2 s )r ξ+1 s sin[(ξ + 1)θ s ] (10) u y (r, θ) = 2K (ξ) II µ √ 2π 2α s α d r ξ+1 d cos[(ξ + 1)θ d ] -(1 + α 2 s )r ξ+1 s cos[(ξ + 1)θ s ] . (11)
The values of ξ used to describe the experimental curves were obtained through the measured evolution of W bd with D as discussed earlier (Brener and Bouchbinder, 2021a). The stress intensity factor was computed as (Eq. ( 5)

from (Brener and Bouchbinder, 2021a)):

K (ξ) II = EW bd (D fin ) (1-ν 2 )f II (C f )r (1+2ξ) fin
, with E, ν respectively the elastic modulus and Poisson's ratio, and

f II (C f ) = αs (1-ν)R(C f ) C 2 f C 2 S
the universal function of rupture velocity.

Description of strain perturbations with theoretical predictions

We now compare the theoretical predictions to experimental strain and displacement evolution of two different frictional ruptures presenting values of ξ = -0.32, -0.27, and final values of W bd of 9.5 and 11 J/m 2 , respectively (Fig. 3b). The values of singularity order were obtained as described in the Results section, by computing the slope of the second branch of the W bd vs D curve and using Eq. 5.

The above-mentioned comparison is presented in Fig. 4. Note that for both models, i.e. LEFM and unconventional theory, the predictions of strain fail ahead of the rupture tip. This is explained by the fact that the two models assume a dynamic rupture driven along an infinite fault by a shear stress equal to the residual stress. As such, they overlook any finite-size effects emerging from the finiteness of the specimen size and the distance to the applied boundary conditions. Moreover, please note that the measurement location was chosen to be the closest possible to the fault plane (strain gauges at ∼1 mm), to capture stress and displacement evolution close to the ones occurring on-fault. However, this choice implies the likelihood of performing measurements within the cohesive zone, expected to be for PMMA around 2-5 mm. The cohesive zone (indicated in Fig. 4a-d with the shaded grey area) was excluded when performing the LEFM fits, given that this model assumes conditions of small-scale yielding (dissipation zone small with respect to the other length scales).

The experimental data were compared with the predictions of LEFM (ξ = -0.5) inverting G c from the best possible fit. The inversion and the minimization algorithm employed to obtain the best solution of G c use simultaneously two strain components (∆ε xx , ∆ε xy ) following the method described in previous studies [START_REF] Svetlizky | Classical shear cracks drive the onset of dry frictional motion[END_REF] (Fig. 4). ∆ε xx and ∆ε xy are obtained by subtracting the initial strain from ε xx and the residual strain from ε xy .

The best fits output values of G c slightly different from the values of W bd estimated through the integration of the slip stress curves. The LEFM predictions do not deviate excessively from the experimental curves for either event, showing an acceptable but not accurate description of the strain perturbations for ∆ε xx and ∆ε xy (Fig. 4a-b). A stronger deviation is observed for ∆ε xy , particularly in the case of ξ = -0.27, independently of the distance from the rupture tip (Fig. 4b). In the second stage, predictions accounting for the unconventional model were computed. The values of ξ and W bd measured as described in the Results section were imposed. The unconventional model returned, for the two events, satisfactory predictions of the evolution of ∆ε xx and ∆ε xy (Fig. 4a,b). It can be noted that the greater the deviation from ξ = -0.5, the greater the disparities between LEFM and the unconventional model (Fig. 4a,b). In addition, the prediction obtained for u x (Fig. 4c,d) is close to the experimental evolution in terms of magnitude. However, while u x evolution is similar within the first microns, the experimental data deviate from the theoretical prediction far behind the rupture tip (Fig. 4c,d). The model returned reasonable predictions of u x for ξ=-0.32, and adequate ones for ξ=-0.27. predictions of the evolution of breakdown work with slip behind the crack tip.

Starting from the stress evolution estimates computed for both LEFM and unconventional model, the breakdown work was computed following Eq. 4.

Please note that neither of the two models is expected to describe the experimental evolution for small D, given that both of them rely on a small-scale yielding condition. Concerning the overall evolution, LEFM predictions deviate in both quantity and temporal evolution from the experimental data. On the contrary, the unconventional model provides a good prediction, particularly for D > D c , as expected from the unconventional theory (Fig. 4e-f).

These results highlight that while LEFM provides reasonable estimates of fracture energy, the unconventional theory provides more coherent predictions of breakdown work evolution with slip when enhanced weakening is observed.

Flash heating as a possible weakening mechanism

These results provide the first complete evidence of unconventional stress fields during the dynamic propagation of laboratory frictional rupture, caused by continuous stress weakening behind the rupture tip. The observed unconventional singularity orders could emerge, among others, from frictional weakening mechanisms such as; thermal activation [START_REF] Bar-Sinai | On the velocity-strengthening behavior of dry friction[END_REF], viscous friction [START_REF] Brener | Frictional shear cracks[END_REF], powder lubrication [START_REF] Reches | Fault weakening and earthquake instability by powder lubrication[END_REF], flash heating [START_REF] Molinari | Dependence of the coefficient of friction on the sliding conditions in the high velocity range[END_REF][START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] Brantut | The fracture energy of ruptures driven by flash heating[END_REF], thermal pressurization [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] Viesca | Ubiquitous weakening of faults due to thermal pressurization[END_REF]. Among these, flash heating has been shown to be activated under similar experimental conditions in Homalite [START_REF] Rubino | Understanding dynamic friction through spontaneously evolving laboratory earthquakes[END_REF], and thus could be the best candidate to explain the unconventional stress fields observed in our experiments. Moreover, the high slip rate measured near-fault promotes the activation of flash heating as previously shown [START_REF] Molinari | Dependence of the coefficient of friction on the sliding conditions in the high velocity range[END_REF][START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] Goldsby | Flash Heating Leads to Low Frictional Earthquake Slip Rates[END_REF]. This agrees with the clear dependence of ξ values with maximum slip rate and rupture velocity observed in our events (Fig. 3c): higher V max are associated with ξ values that deviate from the conventional value (-0.5).

Flash heating is activated when the fault slip velocity becomes higher than a critical weakening slip velocity V w , causing mechanical degradation of contact asperities during their lifetime [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] Goldsby | Flash Heating Leads to Low Frictional Earthquake Slip Rates[END_REF]. The temperature reached at the asperities was computed trough

T asp = T amb + 1 (ρcp √ kπ) τ c V √ t c
with T amb the ambient initial temperature, τ c the stress acting on the single asperity, t c the lifetime of a contact, ρ the bulk density, c p the bulk specific heat and k the thermal diffusivity. Under our experimental conditions, the temperature increased with slip velocity, exceeding the material's melting temperature (T asp > T melting = 160 • ) (Fig. 5a, b), and indicating that melting of asperities probably occurred in our experiments [START_REF] Rubino | Understanding dynamic friction through spontaneously evolving laboratory earthquakes[END_REF]. We compared the evolution of W bd with D, normalized respectively by G c and D c , with asymptotic solutions for flash heating phenomena [START_REF] Brantut | The fracture energy of ruptures driven by flash heating[END_REF].

For D < D c (small slip), the evolution of W bd can be described by the asymptotic solution derived for adiabatic conditions [START_REF] Brantut | The fracture energy of ruptures driven by flash heating[END_REF]:

W bd = ρc(T m -T amb )w √ 2π D V t A w + D 2 ( 12
)
where t A w = ρc(T m -T f )/τ a ( √ 2πw)/V w (time required for a layer of thickness √ 2πw to reach T melting ), w is the shear zone thickness (assumed here as w = 4a with a the asperity size), and τ a is a normal stress dependent contact shear stress at the origin of the change in temperature in the fault layer (Fig. 5c). In presence of gouge along the interface, τ a will correspond to the macroscopic shear stress τ 0 . Along bare rock interfaces,

τ a = τ c a ∆Lasp ,
where ∆L asp is the average distance between two asperities (see Annex A for details). Note that this model assumes a constant sliding velocity V . This assumption looks fairly reasonable in our case, as the first part of the stress weakening (D < D c ) occurs in a very short time window, during which V is nearly constant.

For D > D c , a second asymptotic solution considering the coupled elastodynamics and frictional motions of the propagating rupture can be used [START_REF] Brantut | The fracture energy of ruptures driven by flash heating[END_REF]:

W bd = τ c D SP w µV w 3πτ a C f (1/3) D D SP w (2/3) (13)
where D SP w = V w α( ρc(Tw-T f ) τaVw ) 2 is a characteristic slip weakening distance.

While this asymptotic solution is expected to describe the evolution of breakdown work at a larger seismic slip than the one observed in our experiments, this equation can still be used here because (i) heat diffusion at the scale of asperities is expected to control fault weakening when D > D c and (ii) τ a increases with τ 0 , through the increase of a ∆Lasp with σ n .

Assuming our experimental estimate of C f , this asymptote well describes the second branch of the evolution of W bd with D (power law with an exponent of 2/3, Fig. 5c). Such scaling is also observed at large slip for thermal pressurization in drained conditions, suggesting that this exponent is related to diffusion mechanisms regulating the weakening of faulting during seismic slip [START_REF] Brantut | The fracture energy of ruptures driven by flash heating[END_REF][START_REF] Viesca | Ubiquitous weakening of faults due to thermal pressurization[END_REF].

Importantly, an energy dissipation W bd greater than the fracture energy G c was already observed in [START_REF] Barras | The Emergence of Crack-like Behavior of Frictional Rupture: Edge Singularity and Energy Balance[END_REF] for sliding interfaces whose frictional behavior is described by a rate-and-state friction law. Despite this excess, the rupture dynamics were well described by a conventional LEFM analysis (with ξ = -0.5). This was later justified by Brener and Bouchbinder (2021b), who showed that ruptures along interfaces obeying rate-and-state friction displayed a singularity ξ = -0.406 ≃ -0.5, which corresponds to the lower end of our measurements. However, fault characteristics (e.g. roughness, fluid diffusivity, etc.) and external factors such as initial stress state or on-fault temperature can alter the friction law that controls interface slip (i.e. flash heating, thermal pressurization, and others) and change the singularity observed near the rupture accordingly. In the case of flash heating, the observed evolution of breakdown work with slip generates, for example, a singularity order ξ = -0.25 [START_REF] Brantut | The fracture energy of ruptures driven by flash heating[END_REF], which corresponds to the higher-end exponents of Fig. 3. In our experiments, continuous values of exponents ξ have been measured between ξ = -0.42 (rate-and-state) and ξ = -0.22 (flash heating). This can be caused by the presence of a population of contact asperities, each of which have a different size, experience a different normal and shear stress, and reach thus a different value of temperature and slip velocity during rupture (implying that not necessarily all contact asperities experience flash heating) resulting, on average, in a smooth transition from rate and state frictional contact for the lower slip velocities (nearly conventional, ξ = -0.5) to flash heating for larger slip velocities (unconventional, ξ = -0.25).

Implications and conclusions

These experimental results show that the continuous weakening activated along the fault can modify the singularity order governing displacement and stress fields around the rupture tip, inducing a slip and scale-dependent breakdown work, rather than a constant one. Moreover, this work highlights from an experimental point of view that frictional rupture analysis in the linear elastic fracture framework might not always be sufficient when frictional weakening mechanisms occur away from the rupture tip. Importantly, as long as the residual stress does not reach a steady-state value far from the rupture tip, as happens for thermal weakening processes, the singular fields will hardly recover the conventional square-root singularity, independently of the rupture size. One could nonetheless assess the dynamics of such earthquakes, building on a Griffith criterion adapted to unconventional singularities (see Eq. ( 7) (Brener and Bouchbinder, 2021a)). However, this would involve both the fracture energy and the cohesive zone size that often depends on the structural problem (loading conditions, fault geometry).

Furthermore, the activation of thermal mechanisms depends not only on the rupture characteristics such as crack velocity but also on ambient conditions (such as initial temperature) and possibly slip history controlling asperity roughness and strength. As a result, both rupture dynamics and fault weakening are expected to be governed by fault geometry and rheology and may vary depending on the natural environment.

Our new results highlight the difficulty in a priori estimating the relevant parameters governing the dynamics of the seismic rupture, expected to control the final rupture length (earthquake size). One may legitimately wonder whether theoretical models will be able to capture these complex behaviors, or whether numerical simulations, as proposed in recent studies, will be required instead [START_REF] Lambert | Rupture-dependent breakdown energy in fault models with thermo-hydro-mechanical processes[END_REF].

However, together with the recent development of the unconventional singularity theory (Brener and Bouchbinder, 2021a), our results open the door for a better understanding of the rupture dynamics and energy budget of natural earthquakes, through the possible evaluation of the equations of motions for unconventional rupture phenomena.

Figure 1

 1 Figure 1: a. Experimental setup -Direct shear biaxial apparatus with PMMA samples generating an artificial fault. Strain gauges rosettes are located along the fault at a distance of 1-1.5 mm from the fault plane. b. Temporal evolution of vertical strain (obtained through high-frequency strain gauges acquisition system) at the three different locations along the fault. When the fault experiences instability, the shear rupture propagates along the interface and causes a strain perturbation concurrent with the passage of the front (indicated by the blue arrows). Yellow shaded areas indicate the time window selection shown in the following panel. c. Zoom-in of (b). d. Zoom-in of (.c) The red curve indicates the strain gauge location shown in panel (e). e. Vertical strain temporal evolution for the central location. Please note that the y-axis and x-axis limits change for each panel. t0 refers to the arrival time of the rupture for the central strain gauge.

Figure 2 :

 2 Figure 2: Elastic fields around the rupture tip. Evolution of a. shear stress computed from the measured shear strain ε xy , b. slip velocity computed from the measured horizontal strain ε xx , and c. material displacement computed from the estimated slip velocity for several events presenting different C f (colorbar).

Figure 3 :

 3 Figure 3: Slip-dependent breakdown work and the emergence of unconventional singularities. a. Evolution of (τ -τ res ) with D defining the fault's weakening for different events. The integration of these curves leads to the evolution of W bd with D for different C f (b).

Cs 2 ,

 2 where (C d ,C s ) are respectively the P-wave and S-wave velocity, and R(C f ) = 4α d α s -(1+α 2 s ) 2 is the Rayleigh function. (r, θ) are corrected for the distortion induced by the dynamic rupture velocity C f , becoming θ d = arctan(α d tan(θ)), θ s = arctan(α s tan(θ))

FinallyFigure 4 :

 4 Figure 4: Strain and displacement field described by unconventional singularity for two different events (respectively top and bottom panels). a., b. Comparison of the measured strain perturbations ∆ε xx and ∆ε xy with the theoretical predictions considering: i) the estimated unconventional singularities respectively ξ=-0.32 (a) and ξ=-0.27 (b), and G = W bd (in black) and ii) the LEFM conventional singularity ξ=-0.5 with G = G c (the best fit)(in grey). c., d. Evolution of the material displacement u x with predictions for unconventional and conventional singularity. e., f. Comparison of the experimental evolution of breakdown work with slip estimated at gauge location with theoretical predictions for unconventional theory (black solid line) and LEFM (grey solid line).

Figure 5

 5 Figure 5: a. Evolution of local shear stress τ , with slip velocity for one event. b. Temperature evolution with slip velocity at asperity scale compared with melting temperature of PMMA (T m = 160). c. Slip dependence of breakdown work (curves are normalizedrespectively by G c and D c ). W bd evolution exhibits two power laws with exponents of 2 and 0.6. The experimental curves are all described by the asymptotic solutions related to an adiabatic regime for small D and a diffusive regime for large D[START_REF] Brantut | The fracture energy of ruptures driven by flash heating[END_REF]. The dotted black line shows the expected evolution of W bd assuming LEFM at the strain gauges position.

Appendix A.

For the estimate of ∆L asp , a simplified description of the interface roughness is used, considering only one population of asperities of typical size a and height h, separated by an average distance ∆L asp . The number of asperities was computed considering the following relationship Ar An = Gc


Sergio Jiménez-Manchón

Florent Rivals

Lionel Gourichon

Gabriel De Prado

Ferran Codina

Pere Castanyer

Joaquim Tremoleda

Marta Santos

Armelle Gardeisen

A combined approach to reconstructing livestock management in Iron Age north-eastern Iberia: estimating the season of death and palaeodiet using cementochronology and dental micro-and mesowear analyses

Keywords: PALAEODIET, PALAEOENVIRONMENT, SEASONALITY, EMPÚRIES, ULLASTRET

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A combined approach to reconstructing livestock management in Iron Age north-eastern Iberia: estimating the season of death and palaeodiet using cementochronology and dental micro-and mesowear analyses

Sergio Jiménez-Manchón, Florent Rivals, Lionel Gourichon, Gabriel de Prado, Ferran Codina, Pere Castanyer, Joaquim Tremoleda, Marta Santos, Armelle

Gardeisen

To cite this version:

Sergio Jiménez-Manchón, Florent Rivals, Lionel Gourichon, Gabriel de Prado, Ferran Codina, et al.. A combined approach to reconstructing livestock management in Iron Age north-eastern Iberia: estimating the season of death and palaeodiet using cementochronology and dental micro-and mesowear analyses. Archaeofauna, International Journal of archaeozoology, 2023, 32 (1), 10.15366/archaeo-fauna2023.32.1.010. halshs-04302125

INTRODUCTION

The Greek emporion of Empúries (L'Escala, Girona) was established in the region of Empordà, on the Mediterranean coast of the north-eastern Iberian Peninsula, in the early 6 th century BC (Aquilué et al., 1999;[START_REF] Santos | Emporion arcaica: los ritmos y las fisonomías de los dos establecimientos originarios, a partir de los últimos datos arqueológicos[END_REF]. At Empúries, the economy mainly focused on trade with Iberian populations and other Mediterranean peoples. At the same time, agriculture developed significantly in the surrounding areas (e.g., [START_REF] Francès | Les implicacions econòmiques i comercials de la concentració de reserves de cereals a la Catalunya costanera en època ibèrica[END_REF], generating a series of specific settlements named "silo champs" by archaeologists. The Greek pottery recovered from these sites shows that they were closely related to Empúries. This Greek presence in the Empordà region influenced local populations, as shown by the agglomeration of Ullastret, located 15 km away [START_REF] Santos | Movilidad y contacto de culturas en el extremo nordeste de la península ibérica: lecturas desde la arqueología en el núcleo portuario foceo de Emporion y la ciudad ibérica de Ullastret[END_REF]. The urban planning of this latter city was developed as early as the 6 th century BC. The architecture shows Greek influences and archaeological excavations uncovered a significant number of Greek pottery imports [START_REF] Codina | Les imitations de céramique coloniale des sites ibériques d'Ullastret (Catalogne)[END_REF][START_REF] Santos | Movilidad y contacto de culturas en el extremo nordeste de la península ibérica: lecturas desde la arqueología en el núcleo portuario foceo de Emporion y la ciudad ibérica de Ullastret[END_REF]. All of this makes Ullastret one of the most important indigenous sites in the northeast of the Iberian Peninsula (Martin et al., 2010;[START_REF] Codina | ): balance de una década de investigaciones y nuevas perspectivas[END_REF].

The palaeoenvironment of the Empordà region was characterized by the proliferation of wetlands, predominantly small lagoons and ponds. Palaeonvironmental data highlight a maximum rise in sea level during the Iron Age, leading to the spread of marshy areas and salt pans [START_REF] Julià | Coastal Evolution in a Mediterranean Microtidal Zone: Mid to Late Holocene Natural Dynamics and Human of the Castelló Lagoon, NE Spain[END_REF]. Wooded areas were also widespread, with abun-dant scrubs and holm oaks [START_REF] Montaner | El paleopaisatge fluvio-estuari d'Empúries[END_REF]Castanyer et al., 2016;Ejarque et al., 2022). The reduction of forest mass indicates progressive deforestation during the Iron Age (e.g., [START_REF]Vegetation history and human activity during the last 6000 years on the central Catalan coast (northeastern Iberian Peninsula)[END_REF][START_REF]Paisatge i explotació forestal durant el I mil•leni a. n. e. a la plana empordanesa[END_REF], linked to an increase in arable land [START_REF] Valenzuela-Lamas | Economia i canvi socio-culturala Catalunya durantl'edat del ferro[END_REF].

The presence of fortified and walled cities, as in the case of Ullastret, and the emergence of local elites, suggests increased social differentiation in Iron Age times [START_REF] Sanmartí | Paisatges ibèrics. Tipus d'assentaments i formes d'ocupació del territori a la costa central de Catalunya durant el període ibèric ple[END_REF][START_REF]From Local Groups to Early States: The Development of Complexity in Protohistorican Catalonia[END_REF][START_REF] Sanmartí | Els íbers del nord[END_REF]. This is supported by the presence of warrior equipment, weapons and prestige grave goods at the necropolises of Empúries [START_REF] Vives-Ferrándiz | Las necrópolis de Ampurias, I. Introducción y necrópolis griegas[END_REF] and Puig de Serra (Serra de Daró-Ullastret).

The political, social and economic context had an impact on herding practices. However, herding practices are still not fully understood in the north-eastern Iberian Peninsula (e.g., [START_REF]Noves dades sobre la ramaderia especialitzada en el treball i el transport dins les societats camperoles del bronze final i la primera edat del ferro a les depressions de l'Empordà i la Selva (Girona)[END_REF][START_REF] Valenzuela-Lamas | Shipping amphorae and shipping sheep? Livestock mobility in the north-east Iberian peninsula during the Iron Age based on strontium isotopic analyses of sheep and goat tooth enamel[END_REF]. Territorialisation and an increase in cultivated areas probably led to a reduction in grazing areas and, consequently, to risks of animal malnutrition, as suggested by [START_REF] Valenzuela-Lamas | Animal Husbandry Across the Western Roman Empire. Changes and Continuities[END_REF]. Several authors have postulated that agriculture and animal husbandry were in competition with each other [START_REF] Valenzuela-Lamas | Alimentació i ramaderia al Penedès durant la protohistòria (segles VIII-III aC). Premi d'Arqueologia[END_REF][START_REF]La gestió dels animals al nordest de la Península Ibèrica entre els segles V ANE. -V DNE. Proposta Metodològica d'Integració de les anàlisis arqueozoològiques als estudis de cronologies històriques[END_REF]. Colominas et al. (2011) proposed the hypothesis of differential management of the environment and grazing areas according to the stockbreeding species in the Empordà region. More specifically, they suggest that Caprinae (sheep and goats), which have a flexible diet, were kept in wooded areas or in scrublands. Cattle generally need more water and fresh pastu-diferentes en la misma región -y en un mismo entorno-hace que nos preguntemos por la variabilidad de las prácticas ganaderas en estos asentamientos, así como su relación con el entorno. Para ello, además del análisis zooarqueológico tradicional, hemos reconstruido la paleodieta de los animales y la estacionalidad de muerte. Para ello, se han combinado tres proxies: el análisis del meso-y microdesgaste dental, y la cementocronología. Los análisis realizados en dientes de caprinos y bovinos de los yacimientos de Empúries y Ullastret (s. VI-IV a.C.) apoyan la hipótesis propuesta en trabajos anteriores sobre el uso diferente del medio según las especies: los caprinos (ovejas y cabras) se alimentaron en áreas marginales y menos aptas para la agricultura. Por el contrario, el ganado vacuno se alimentó en pastos y en zonas húmedas. La combinación entre la cementocronología y el microdesgaste sugiere un consumo menor de herbáceas en los meses de verano. Y en cuanto a la estacionalidad, los resultados han mostrado que los caprinos fueron principalmente sacrificados entre los meses de primavera y verano.

PALABRAS CLAVE: PALEODIETA, PALEOAMBIENTE, ESTACIONALIDAD, EMPÚRIES, ULLASTRET, NORESTE DE LA PENÍNSULA IBÉRICA res than Caprinae and would thus have been kept in grassland and wetland areas. This hypothesis has also been put forward for the neighbouring region of Languedoc in France (Nieto-Espinet et al., 2020a).

The reconstruction of the animal diet can shed some light on this debate. However, studies of the animal diet are scarce for such areas (e.g., [START_REF] Gallego-Valle | Pastando en las marismas. Una aproximación desde la técnica del desgaste dentario a la alimentación del ganado ovino en el Empordà romano (noreste de la Península Ibérica)[END_REF][START_REF] Jiménez-Manchón | Discriminating management strategies in modern and archaeological domestic caprines using low-magnification and confocal dental microwear analyses[END_REF]. Here, we aim to provide new data on the palaeodiet of livestock species from Ullastret and Empúries. To this end, we used two methods to obtain relevant economic as well as environmental information at different time scales: dental mesowear (i.e., the diet of the last years of the animal's life) and dental microwear (i.e., diet of the last weeks of the animal's life) analyses. These proxies can also distinguish different intra-and inter-species feeding strategies. The seasonality of pastoral practices is a key element for reconstructing herding practices, but is still poorly known. Cementochronology can potentially shed light on this grey area since it provides information on the season of death through the study of incremental growth marks in the cementum. In addition, the combination of dental mi-crowear analysis and cementochronology can characterize the type of diets in relation to the season of death. Here, we studied the cementochronology of livestock species from both sites.

In this work, we present the (1) dental microwear analysis, (2) dental mesowear analysis, (3) cementum analysis, and (4) zooarchaeological analysis (i.e., taxonomical representation and mortality profiles) of unpublished assemblages from Ullastret and Empúries. This multiproxy approach aims to provide new data to enhance our understanding of the economic role of animal husbandry and landscape use in Iron Age Empordà.

MATERIALS AND METHODS

The archaeological sites

The faunal samples come from two Iron Age sites in the Empordà region in Spain: Ullastret (Girona) and Empúries (L'Escala, Girona) (Figure 1). Empúries (from the Greek emporion) consisted FIGURE 1 Location of Empúries and Ullastret in the northeast of the Iberian Peninsula.

Archaeofauna 32-1 (2023): 00-00 of two settlements: Palaiapolis ('the old town' in Greek) founded around 580 BC on a semi-island, and Neapolis ('the new town' in Greek) founded around 525-500 BC (Aquilué et al., 1999;[START_REF] Santos | Emporion arcaica: los ritmos y las fisonomías de los dos establecimientos originarios, a partir de los últimos datos arqueológicos[END_REF][START_REF] Castanyer | Nuevos datos arqueológicos sobre la evolución urbana de Emporion[END_REF]. The studied faunal remains derive from phases IIId and IIIe of Palaiapolis, henceforth Empúries B, dated to the early Iron Age I (550-450 BC), and from habitation B in square 21 at Neapolis. Faunal remains from the late Iron Age I (450-325 BC), henceforth Empúries C, come from two assemblages from Neapolis: Stratigraphic Units (SU) 5016 and 6204. SU 5016 (400-375 BC) is located in the northern part of the Temenos and yielded a significant number of culinary pottery sherds [START_REF]Datación de la muralla griega meridional de Ampurias y caracterización de la facies cerámica de la ciudad en la primera mitad del siglo IV a[END_REF][START_REF] Delgado | ¿Dualidad étnica o heterogeneidad social? Equipos cerámicos y prácticas cotidianas en la Neápolis de Empotion, c. 425-375 a[END_REF]. SU 6204 is a large deposit of pottery and faunal remains located in the northern sector of the temple of Asklepeion, dedicated to Asklepios, located in the south of the city.

The archaeological site of Ullastret is an Iberian city located just 15 km south of Empúries. The site consists of two fortified complexes separated by a palaeo-pond: Puig de Sant Andreu, located on a small hill (50 m. a.s.l.), and Illa d'en Reixac, located on a plain and surrounded by the same palaeo-pond. This 'binary' settlement, with a total area of almost 15 ha during the 4 th century B.C, is the most complex indigenous settlement in the northeast of the Iberian Peninsula (Martin et al., 1999(Martin et al., , 2010)). Faunal remains dated to the late Iron Age I (550-450 BC, phase Ullastret II/III) were found in the levels of insula 7 of Illa d'en Reixac (Martin et al., 1999). The early Iron Age II material (phase Ullastret IV/V phase) comes from the insulae 7 and 5 of Illa d'en Reixac.

Zooarchaeological study

The reference collection housed at the ASM/ UMR5140 laboratory (CNRS, Montpellier) and various atlases and catalogues of comparison [START_REF] Lambert | Atlas ostéologiques pour servir à l'identification des mammifères du Quaternaire. II: herbívores: tête, rachis, ceintures scapulaire et pelvienne, membres[END_REF][START_REF] Schmid | Atlas of Animal Bones[END_REF][START_REF] Barone | Anatomie comparée des mammifères domestiques[END_REF][START_REF]Diferencias osteológicas entre las ovejas (Ovis aries Linné) y cabras (Capra hircus Linné)[END_REF][START_REF] Helmer | Discrimination des genres Ovis et Capra à l'aide des prémolaires inférieures 3 et 4 et interprétation des âges d'abattage: l'exemple de Dikili Tash[END_REF][START_REF] Halstead | Sorting the sheep from the goats: morphological distinction between the mandibles and mandibular teeth of adult Ovis and Capra[END_REF][START_REF] Wilkens | Archeozoologia[END_REF][START_REF]Assessing the reliability of criteria used to identify mandibles and mandibular teeth in sheep, Ovis, and goats, Capra[END_REF][START_REF]Assessing the reliability of criteria to identify postcranial bones in sheep, Ovis, and goats, Capra[END_REF] were used for the anatomical and taxonomical identification of the remains. When the taxon could not be determined, bones were classified according to body mass categories, taking into account the faunal spectrum of the assemblage and age. In this work, 'large-sized' mammals (<150 kg) include adult cattle (Bos taurus), red deer (Cervus elaphus) and equids (Equus sp.); 'medium-sized' (15-100 kg) include immature cattle, equids and red deer, as well as adult sheep (Ovis aries), goats (Capra hircus), suid (Sus sp.) and dog (Canis familiaris); and finally, 'small-sized' animals (>15 kg) comprise immature sheep, goat, pig, dog, small carnivores and lagomorphs. Age-at-death was estimated on the basis of the degree of epiphysation and tooth eruption and wear [START_REF] Payne | Kill-off paterns in sheep and goats: The mandibles from Asvankale[END_REF][START_REF] Grant | The use of the tooth wear as a guide to the age of domestic animals[END_REF][START_REF]Exploitation des prélèvements et fichiers de spécialité (PRL, FAUNE, OS)[END_REF]. The number of identified specimens (NISP), the minimum number of elements (MNE) and the minimum number of individuals (MNI) [START_REF] Lyman | Vertebrate taphonomy[END_REF] were calculated to assess the integrity of the assemblage. Mortality profiles were constructed using the Bayesian method established by [START_REF] Valenzuela-Lamas | Assessment and prediction of mortality profiles. A Bayesian approach -Proposta i predicció de models d'aprofitament ramader. El cas de la Cossetània oriental entre els segles VII-III aC[END_REF] and interpreted according to the zootechnical models of animal production strategies described by several authors [START_REF] Stein | Regional Economic Integration in Early State Societies. Third Millenium B.C. Pastoral Production at Gritille, Southeast Turkey[END_REF][START_REF] Helmer | La gestion des caprinés domestiques dans le Midi de la France[END_REF][START_REF] Blaise | L'élevage au Néolithique dans le sud-est de la France : éléments de réflexion sur la gestion des troupeaux[END_REF][START_REF] Helmer | 00-00 sudest de la France : saisonnalité des abattages,rélations entre grottesbergeries et sites de plein air[END_REF][START_REF] Helmer | The development of the exploitation of products from Capra and Ovis (meat, milk and flece) from the PPNB to the Early Bronze in the northern Near East (8700 to 2000 BC cal.)[END_REF]. A total of 3525 faunal remains from Empúries (N= 2012) and Ullastret (N= 1513) were analysed (Table 1). Number of specimens analysed by archaeological and site according to the approach in this study.

Dental microwear analysis

Dental microwear analysis was carried out following the method defined by Solounias & [START_REF] Semprebon | Advances in the reconstruction of ungulate eco morphology with application to early fossil equids[END_REF]. The occlusal surface of teeth was carefully cleaned with acetone and then with 96% ethanol [START_REF] Godfrey | Can low-magnification stereomicroscopy reveal diet?[END_REF]. Afterwards, the occlusal surface was moulded using the high-precision polyvinylsiloxane dental impression material Provil Novo Light C.D.2 (Heraeus Kulzer) for the first contact layer, and Putty Blue Eco (Detax) for the supporting layer. These moulds were used to make transparent epoxy resin casts (C.P. Química CPOX P 1069/A -CPEN 1585/B). The epoxy casts were analysed using a stereomicroscope at low magnification (x35). For each tooth, we analysed two different areas of 0.4 mm x 0.4 mm Archaeofauna 32-1 (2023): 00-00 (0.16 mm 2 ) and averaged the data on the protoconid and/or hypoconid of the first (M1), second (M2) or third (M3) lower molars, and on the paracone and/or metacone of the upper M1, M2 and/or M3 [START_REF] Xafis | Which tooth to sample? A methodological utility of premolar/ non-carnassial teeth in microwear analysis of mammals Palaeogeography[END_REF]. Analyses were performed on well-preserved enamel surfaces, avoiding taphonomic alterations such as fissures, damaged areas and trampling features [START_REF] King | Effect of Taphonomic Processes on Dental Microwear[END_REF][START_REF] El-Zaatari | Occlusal microwear texture analysis and the diets of historical/prehistoric hunter-gatherers[END_REF]Uzunidis et al., 2021). Five microwear variables were recorded: the number of pits (NP), the number of scratches (NS), the presence/absence of gouges (G), the scratch width score (SWS) and the presence/absence of cross scratches (XS). We then generated a bivariate plot graph featuring the total number of scratches and pits using the R Statistical Software (version 3.5.1) and the R code from [START_REF]Dietary traits of ungulates in northeastern Iberian Peninsula: Did these Neanderthal preys show adaptive behaviour to local habitats during the Middle Palaeolithic?[END_REF]. After discarding teeth with taphonomic alterations [START_REF] King | Effect of Taphonomic Processes on Dental Microwear[END_REF][START_REF] Rivals | Effect of Ontogenetic-age Distribution in Fossil and Modern Samples on the Interpretation of Ungulate Paleodiets Using the Mesowear Method[END_REF][START_REF] El-Zaatari | Occlusal microwear texture analysis and the diets of historical/prehistoric hunter-gatherers[END_REF]Uzunidis et al., 2021), we analysed a total of 131 caprine and cattle teeth from Empúries (N= 58) and Ullastret (N= 73, Table 1). The teeth sampled for dental microwear come from the faunal remains studied in this work and from other faunal assemblages analysed by other researchers (Aquilué et al., 1999;[START_REF]Noves dades sobre la ramaderia especialitzada en el treball i el transport dins les societats camperoles del bronze final i la primera edat del ferro a les depressions de l'Empordà i la Selva (Girona)[END_REF].

Dental mesowear analysis

The Mesowear Score (MWS) consists in the macroscopic and qualitative analysis of the morphology of the dental cusps located on the buccal side of upper molars and on the lingual side of lower molars [START_REF] Fortelius | Functional characterization of ungulate molars using the abrasion-attrition wear gradient: a new method for reconstructing paleodiets[END_REF][START_REF] Franz-Odendaal | Differential mesowear in the maxillary and mandibular cheek dentition of some ruminants (Artiodactyla)[END_REF][START_REF] Kaiser | Extending the tooth mesowear method to extinct and extant equids[END_REF]. Cusp morphology is classified into seven categories, where '0' is characterized by high relief and pointed cusps, corresponding to the least abrasive wear, and '6' is characterized by low relief and flattened (blunt) cusps, corresponding to the most abrasive wear [START_REF] Rivals | Effect of Ontogenetic-age Distribution in Fossil and Modern Samples on the Interpretation of Ungulate Paleodiets Using the Mesowear Method[END_REF]. A browser dietary regime varies between a MWS value of 0 to 2.00, a grazer regime between 2.09-5.47, and mixed feeders between 0.40-2.74 [START_REF] Fortelius | Functional characterization of ungulate molars using the abrasion-attrition wear gradient: a new method for reconstructing paleodiets[END_REF][START_REF] Rivals | Effect of Ontogenetic-age Distribution in Fossil and Modern Samples on the Interpretation of Ungulate Paleodiets Using the Mesowear Method[END_REF][START_REF]Dietary change and evolution of horses in North America[END_REF]. A total of 104 caprine and cattle teeth from Empúries (N= 56) and Ullastret (N= 48) were included in the analysis (Table 1). The sampled teeth come from the faunal assemblages analysed in this work and from other faunal assemblages studied by other zooarchaelogists (Martin et al., 1999).

Cementum analysis

Dental cementum analysis focuses on the study of incremental growth marks in the cementum of mammalian teeth in order to estimate the age and season at death of the animal [START_REF]Seasonality studies[END_REF][START_REF]Recording structures of Mammals[END_REF][START_REF] Janin | Faune et saisonnalité : l'organisation temporelle des activités de subsistance dans l'Epipaléolithique et le Néolithique précéramique du Levant Nord (Syrie)[END_REF][START_REF] Rendu | La cémentochronologie[END_REF]Naji et al., , 2022)). The dual characteristics of cementum are that it develops throughout the animal's life, i.e., from the formation of the root until tooth shedding, and its deposition rate follows an annual cycle (cold and warm seasons) [START_REF] Muñoz-Cobo | Seasonal deposition of cementum in first lower molars from Cervus elaphus hispanicus[END_REF]. Cementum is a bone tissue deposited on the exterior part of the tooth root, although it can also extend to the enamel (coronal cementum) in hypsodont species [START_REF] Yamamoto | Historical review of the human celular cementum with special reference to an alternating lamellar pattern[END_REF]. The primary function of this tissue is to anchor the tooth root to the alveolar bone via periodontal ligaments to hold the crown during occlusion and chewing processes. This research focuses on the acellular extrinsic fibre cementum (henceforth, AC). The AC is characterized by a successive and regular alternation of parallel bands. The biological mechanisms underlying this alternation of cementum increments are still not fully understood, but several studies have shown that they are correlated with seasonal changes and other biochronological processes [START_REF]Age Determination of Mammals from Annual Layers in Teeth and Bones[END_REF]Grue & Jensen, 1979;[START_REF] Lieberman | The biology of cementum increments (with an archaeological application)[END_REF]Naji et al., 2022). Two types of AC can be distinguished: 'growth zones' (GZ) and 'annuli'. GZ are characterized by fast growth and a low mineralization rate and build up during the warm season. 'Annuli' are characterized by slow growth and a high mineralization rate and form during the cold season [START_REF] Lieberman | The Biological Basis for Seasonal Increments in Dental Cementum and their Application to Archaeological Research[END_REF][START_REF] Janin | Faune et saisonnalité : l'organisation temporelle des activités de subsistance dans l'Epipaléolithique et le Néolithique précéramique du Levant Nord (Syrie)[END_REF][START_REF] Rendu | La cémentochronologie[END_REF]Naji et al., , 2022)). These bands can be microscopically differentiated on the basis of differences in mineralization rates and sometimes in the orientation of the Sharpey fibres [START_REF]The use of cement layers for absolute ageing of mammalian teeth: a selective review of the literature, with suggestions for further studies and alternative applications[END_REF][START_REF] Lieberman | The Biological Basis for Seasonal Increments in Dental Cementum and their Application to Archaeological Research[END_REF][START_REF] Rendu | Planification des activités de subsistance au sein du territoire des derniers Moustériens. Cémentochronologie et approche archéozoologique de gisements du Paléolithique moyen (Pech de-l'Azé I[END_REF].

Season at death can be estimated by determining whether the last cementum band is an 'annulus' (cold season) or a GZ (warm season) [START_REF] Burke | Histological observations of cementum growth in horse teeth and their application to archaeology[END_REF]. The GZ can provide further detail by comparing the relative growth volume of the last deposit to the thickness of the previous complete bands: 'beginning' (0-33.3%), 'middle' (33.3-66.6%) or 'last' (66.6-100%) part of a GZ [START_REF] Janin | Faune et saisonnalité : l'organisation temporelle des activités de subsistance dans l'Epipaléolithique et le Néolithique précéramique du Levant Nord (Syrie)[END_REF][START_REF] Rendu | Planification des activités de subsistance au sein du territoire des derniers Moustériens. Cémentochronologie et approche archéozoologique de gisements du Paléolithique moyen (Pech de-l'Azé I[END_REF]Sánchez-Hernández et al., 2020). The estimation can be approximately related to the respective following seasons: spring, summer and fall [START_REF]Recording structures of Mammals[END_REF]. 1). The teeth were first extracted from the alveolar bone. Then, the teeth, except the crown, were placed in plastic boxes and covered with transparent epoxy resin (DBF epoxy). A total of three thin sections were cut per root using a slow-speed diamond saw (Buehler IsoMet 1000) in the CEPAM laboratory (CNRS, Nice). After cleaning, the sections were glued (Geofix epoxy) to a glass plate (25x30 mm). The inner face of these sections was mechanically abraded to a thickness of 50-80 μm using a diamond grinding wheel (Buehler PetroThin) and a thin glass coverslip was then glued directly onto the surface of the sample. Thin sections were analysed using a polarizing light microscope at x100 and x200 magnifications. Teeth with taphonomic alterations were discarded [START_REF] Rendu | Approche taphonomique en Cémentochronologie : réexamen du niveau 4 du Pech-de-l'Azé I (Carsac, Dordogne, France)[END_REF][START_REF] Janin | Faune et saisonnalité : l'organisation temporelle des activités de subsistance dans l'Epipaléolithique et le Néolithique précéramique du Levant Nord (Syrie)[END_REF]. The incremental bands, and particularly the last one, were identified using three optical filters: plane-polarized light, cross-polarized light and full-wave retardation plate (λ plate). Microphotographs of the main regions of interest were systemically taken for archiving and further study.

RESULTS

Zooarchaeological data

A total of 1512 faunal remains from Ullastret II/III (550-450 BC) and Ullastret IV/V (450-325 BC) were studied. For the Ullastret II/III phase, 299 remains were analysed, of which 132 (44.1%) were anatomically and taxonomically identified. The most frequent species are Caprinae (Ovis aries and Capra hircus), followed by pigs (Sus domesticus) and cattle (Bos taurus) (Figure 2A). To Archaeofauna 32-1 (2023): 00-00 a lesser extent, dog (Canis familiaris) and rabbit (Oryctolagus cuniculus) remains were also identified. For the Ullastret IV/V phase, 1213 faunal remains were examined, of which 838 (69%) were identified. Caprinae are the most frequent taxon, followed by pigs and cattle (Figure 2B). Dog, red deer (Cervus elaphus), rabbit and an unidentified carnivore were also recovered. The mortality profiles of Caprinae from both phases show that they were mostly slaughtered before the age of 24 months, followed by animals 24-48 months old. According to the zootechnical models currently used in zooarchaeology (cf. section 2.2), this pattern corresponds to flock management strategies focused on meat production. The mortality profile of cattle from Ullastret II/III suggests that slaughtering centred on animals aged 24-48 months, i.e., in the prime of their age and at the meat/weight optimum. For the second phase (Ullastret IV/V), cattle were mainly slaughtered at 12-24 months old for meat exploitation, and at 48 months old or more in the case of dairy exploitation and the use of animals for traction/ploughing. A total of 2012 faunal remains from the Empúries site were studied. For the Empúries B phase, 500 remains out of 1101 (45.4%) were anatomically and taxonomically identified. The most abundant taxa are Caprinae, followed by cattle and pigs (Figure 2C), and to a lesser extent, dog and rabbit remains. In the Empúries C phase, Caprinae, pig and cattle are the most frequent taxa. Other identified taxa are equids (Equus sp.), dog, red deer and rabbit (Figure 2D). The kill-off profiles of Caprinae in both phases (Empúries B and C) indicate a focus on meat production (6-12 and 12-24 months). The slaughtering of individuals at 24-48 and 48-72 months old also reflects milk exploitation. Two culling peaks were observed for cattle at 12-24 and 24-48 months old in Empúries B, indicating mainly meat production. Conversely, in Empúries C, the main peak is observed at 72-96 months old, which may possibly be related to milk production and traction.

Dental mesowear score

The overall results obtained from the dental mesowear analysis are presented in Figure 3a and Table 2. For Ullastret II/III, sheep and goats were grouped together. In the remaining assemblages, sheep and goats were distinguished. The mesowear score (MWS) of Caprinae from both phases (550-450 BC and 450-325 BC) and sites varies between 1.11 and 2.12. These values indicate a dietary pattern compatible with mixed feeders with a browsing tendency. Only the sheep from Empúries B show values compatible with grazers (MWS= 2.69). Concerning variability between sheep and goats, sheep present a higher MWS than goats. The degree of abrasion is higher for cattle than for Caprinae, varying between 2 and 2.5, indicating a grazing dietary regime. In general, these data suggest the presence of more abrasive particles in resources eaten by cattle than those consumed by Caprinae. Caprine MWS values show higher abrasion in the second phase (450-325 BC) than in the first (550-450 BC), suggesting the consumption of less shrubby plants in the first phase.

Dental microwear pattern

A total of 131 goat, sheep and cattle teeth were analysed after discarding specimens with taphonomically altered enamel surfaces.

For the first chronological phase (550-450 BC) of Ullastret and Empúries, the microwear pattern of sheep and goats presents an intermediate number of pits (between 24 and 28) and an intermediate to high number of scratches (16.0-20.5). These data indicate a mixed feeding dietary regime. The average number of scratches is slightly higher in sheep than in goats in both sites ( Archaeofauna 32-1 ( 2023): 00-00 Archaeofauna 32-1 ( 2023): 00-00 of pits (15.5-19.2) and a slightly higher number of scratches (18-21.8) than Caprinae.

For the second period (450-325 BC), Caprinae present an intermediate number of pits (23.0-26.5) and an intermediate to high number of scratches (18)(19)(20). This dental microwear pattern is compatible with a mixed dietary regime. Caprinae from Empúries show a higher number of scratches than those from Ullastret. Cattle teeth bear a higher number of scratches (between 21 and 22) and a lower number of pits (19)(20)(21)(22) than Caprinae.

Seasonal pattern

A total of 14 caprine and five cattle teeth were analysed. Mixed cementum was found in the upper part of the roots of all teeth, and it was thus difficult to identify and study optimal regions of acellular cementum. In this regard, the presence of cementocytes throughout the cementum layer was observed in 14 specimens. For three caprine teeth from the first phase of Ullastret (550-450 BC), the last cementum deposit shows the beginning of a GZ, suggesting slaughtering in spring. In another caprine tooth, a half-thick GZ (ca. summer) was identified. For Empúries, the nature of the last cementum increment is variable: one individual with the beginning of a GZ (ca. spring), a second with a half-thick GZ (ca. summer), and another with an almost-complete GZ (ca. fall). In the case of cattle (Figure 4), one specimen from Ullastret shows the beginning of a GZ, one from Ullastret a middle GZ, and another a final annulus. For the second period (450-325 BC), the cementum analysis of caprine teeth from Ullastret shows slaughtering in the middle (N=2) and at the beginning (N=1) of the warm season. Observa-tions at Empúries are also variable: two Caprinae were killed in summer and two others during the cold season (winter). The two analysed cattle specimens from Ullastret indicate that one was slaughtered in the middle of the warm season (ca. summer) and the second at the end of this period (ca. fall).

Combination of dental cementum and microwear data

it was possible to combine the data from dental microwear and cementum analyses for 14 teeth: five sheep, seven goats and two cattle. The results, presented in Table 4, show a lower number of scratches in summer (middle of the warm season) compared to the other seasons. Summary of the combination between dental microwear and cementum analysis. The average of the number of scratches according the season of death is indicated (N= number of specimens).

DISCUSSION

The zooarchaeological analysis showed that Caprinae were the most abundant animals in the two studied settlements, in keeping with trends observed in other sites from the same chronological and geographical context (e.g., [START_REF]La gestió dels animals al nordest de la Península Ibèrica entre els segles V ANE. -V DNE. Proposta Metodològica d'Integració de les anàlisis arqueozoològiques als estudis de cronologies històriques[END_REF][START_REF]La gestió dels animals al nordest de la Península Ibèrica entre els segles V ANE. -V DNE. Proposta Metodològica d'Integració de les anàlisis arqueozoològiques als estudis de cronologies històriques[END_REF][START_REF] Nieto-Espinet | Canvis en l'alimentació càrnia a Catalunya entre els segles XII i II aC[END_REF][START_REF] Valenzuela-Lamas | Economia i canvi socio-culturala Catalunya durantl'edat del ferro[END_REF]. At Ullastret, pigs are well represented in both periods, Ullastret IV/V (5 th -4 th c. BC) 32 19.2 5.8 21.8 4.7 0.9 60 0 15 28.9 8 18.5 3.5 0.7 10 0 14 24.1 8.9 16.7 5.7 0.4 10 0

Empúries B (6 th -5 th c. BC) Archaeofauna 32-1 (2023): 00-00 Archaeofauna 32-1 (2023): 00-00 while at Empúries, an increase in the frequency of pigs was observed during the Iron Age II phase (450-325 BC) compared to the previous period (550-450 BC). This trend was also observed in the north-eastern Iberian Peninsula (e.g., [START_REF] Valenzuela-Lamas | Alimentació i ramaderia al Penedès durant la protohistòria (segles VIII-III aC). Premi d'Arqueologia[END_REF][START_REF]La gestió dels animals al nordest de la Península Ibèrica entre els segles V ANE. -V DNE. Proposta Metodològica d'Integració de les anàlisis arqueozoològiques als estudis de cronologies històriques[END_REF]Nieto-Espinet, 2012, 2020a, 2020b;[START_REF] Belarte | Animal husbandry and hunting practices in Hispania Tarraconensis: An overview[END_REF], and in neighbouring regions [START_REF]Exploitation des prélèvements et fichiers de spécialité (PRL, FAUNE, OS)[END_REF][START_REF] Belarte | Animal husbandry and hunting practices in Hispania Tarraconensis: An overview[END_REF]. The increase in arable land in that time (e.g., [START_REF] Francès | Les implicacions econòmiques i comercials de la concentració de reserves de cereals a la Catalunya costanera en època ibèrica[END_REF] could explain the increase in pigs, as they do not require grazing lands and are not in competition with agriculture. Higher frequencies of this taxon could also be explained by an increase in urbanisation and intensive meat production [START_REF] Valenzuela-Lamas | Alimentació i ramaderia al Penedès durant la protohistòria (segles VIII-III aC). Premi d'Arqueologia[END_REF][START_REF]La gestió dels animals al nordest de la Península Ibèrica entre els segles V ANE. -V DNE. Proposta Metodològica d'Integració de les anàlisis arqueozoològiques als estudis de cronologies històriques[END_REF]Nieto-Espinet et al., 2020a).

The mortality profiles of Caprinae reflect a mixed production of milk and meat at the two studied sites and periods. A focus on adult individuals was observed in other Iron Age sites in Catalonia and Languedoc (Valenzuela-Lamas, 2008;Albesso et al., 2013;Nieto-Espinet et al., 2020a;[START_REF] Jiménez-Manchón | Pratiques d'élevage entre l'Empordà et le Languedoc à l'âge du Fer. Archéozoologie[END_REF][START_REF] Jiménez-Manchón | Pratiques d'élevage entre l'Empordà et le Languedoc à l'âge du Fer. Archéozoologie[END_REF]. The slaughtering of adult specimens could be linked to the exploitation of secondary products, such as wool and milk. In this respect, archaeological evidence, such as the finding of pondera (measure of weight), in these sites demonstrates the use of wool and related activities [START_REF] Castro Cruel | Piezas toroides de arcilla en yacimientos ibéricos[END_REF]. However, mortality profiles do not follow a single trend, as they vary depending on settlement types and the chronological period. Thus, at Mas Castellar de Pontós, located in the Empordà region [START_REF]La gestió dels animals al nordest de la Península Ibèrica entre els segles V ANE. -V DNE. Proposta Metodològica d'Integració de les anàlisis arqueozoològiques als estudis de cronologies històriques[END_REF], or La Ramasse in Languedoc (Jiménez-Manchón, 2020), adults are predominant, while in Lattara, young specimens are more frequent [START_REF] Jiménez-Manchón | Pratiques d'élevage entre l'Empordà et le Languedoc à l'âge du Fer. Archéozoologie[END_REF]Nieto-Espinet et al., 2020a, 2020b).

Higher variability is observed among cattle. In the first phase (550-450 BC), milk and meat production are balanced in Empúries, as at the same period in Lattara (Languedoc, France) (Nieto-Espinet et al., 2020a). In Ullastret, slaughtering focused on fully-grown animals (24-48 months), as in la Monédière (Languedoc, France) (Nieto-Espinet et al., 2020a). In the second period (450-325 BC), the mortality profiles of both studied sites show an increase in adult individuals, suggesting a greater exploitation of secondary products, such as milk, or the use of animals for traction/ploughing. These data suggest the use of cattle for dairy products, in addition to Caprinae, to the detriment of high-qua-lity meat. This trend has also been observed in other Iron Age sites in Catalonia and Languedoc (e.g., Valenzuela-Lamas, 2008;Nieto-Espinet et al., 2020a). Pottery sherds found in Iron Age sites in the neighbouring region of Languedoc may possibly be related to milk production [START_REF]Les Gaulois du Midi, de la fin de l'Âge du Bronze à la conquête romaine[END_REF]. Pigs were slaughtered to obtain tender meat (12-18 months) and maximum meat/fat yields (18-36 months). A higher number of animals older than 36 months were only observed in the first phase of Ullastret (550-450 BC).

The dental micro-and mesowear analyses results revealed the overall dietary patterns of sheep, goat and cattle at the Empúries and Ullastret sites during the Iron Age. In general, they show that sheep and goats were mixed feeders with a tendency to browse. This suggests a diet rich in eudicotyledonous plants, shrubs and bushes, with a reduced consumption of graminaceous plants. This type of diet has been recorded in sheep and goats in the Bronze Age Balearic Islands [START_REF] Valenzuela-Suau | Caprine dental microwear reveals livestock management and exploitation of landscape during the Middle and Late Bronze Age of the Balearic Islands (ca. 1500-850 cal. BC)[END_REF], and in the Iberian Peninsula during the Roman [START_REF] Gallego-Valle | Pastando en las marismas. Una aproximación desde la técnica del desgaste dentario a la alimentación del ganado ovino en el Empordà romano (noreste de la Península Ibérica)[END_REF], Iron Age [START_REF] Jiménez-Manchón | Reconstruction of caprine management and landscape use through dental microwear analysis : the case of the Iron age site of el Turó de la Font de la Canya (Barcelona, Spain)[END_REF] and Neolithic periods [START_REF] Sierra | Sheep husbandry in the early Neolithic of the Pyrenees: New data on feeding and reproduction in the cave of Chaves[END_REF]. The diet of sheep from Ullastret and Empúries was slightly more abrasive than that of goats, perhaps due to a sheep preference for herbaceous plants (e.g., [START_REF] Treziny | Comparison and overlap of sympatric wild ungulate diet in Cazorla, Segura and Las Villas Natural Park[END_REF]. This trend has been identified in various archaeological contexts using stable isotopes (e.g., [START_REF]Reconstructing dietary and environmental history from enamel isotopic analysis: time resolution of intra-tooth sequential sampling[END_REF] and dental microwear analyses [START_REF] Jiménez-Manchón | Discriminating management strategies in modern and archaeological domestic caprines using low-magnification and confocal dental microwear analyses[END_REF]. Unlike Caprinae, cattle present a grazing dietary regime, compatible with a diet rich in grasses. These data suggest pasturing on meadows and/ or in less wooded areas. The available palaeoenvironmental data describe a landscape dominated by grasslands with the existence of wooded areas. The results obtained from dental meso-and microwear analyses support the hypothesis put forward by Colominas et al. (2011) regarding landscape use and husbandry practices. They postulate that Caprinae, characterized by a more flexible diet than cattle, would have fed in wooded areas less suitable for agriculture. In contrast, cattle require more fresh pastures than Caprinae, and probably fed in grassland areas. However, the consumption of grasses in the form of fodder by cattle cannot be ruled out, although this type of diet is difficult to identify with dental microwear analyses [START_REF] Gallego-Valle | What is on the menu today? Creating a microwear reference collection through a controlled-food trial to study feeding management systems of ancient agropastoral societies[END_REF]. No evidence of animal pens, where animals can be fodder fed, has been found in Ullastret and Empúries. However, such practices have been Archaeofauna 32-1 (2023): 00-00 recognized in other Iron Age sites from Catalonia [START_REF] Jiménez-Manchón | Pratiques d'élevage entre l'Empordà et le Languedoc à l'âge du Fer. Archéozoologie[END_REF] and cannot be excluded.

The results obtained from the dental mesowear analysis show that Caprinae and cattle consumed less woody plants in the second phase (450-325 BC) than in the first phase (550-450 BC). This scenario is consistent with the progressive reduction of the forest mass documented by palaeoenvironmental studies in the north-eastern Iberian Peninsula (e.g., [START_REF]Vegetation history and human activity during the last 6000 years on the central Catalan coast (northeastern Iberian Peninsula)[END_REF][START_REF]Paisatge i explotació forestal durant el I mil•leni a. n. e. a la plana empordanesa[END_REF].

Using a combination of dental microwear and cementum analyses, we characterized the type of animal diet and the season of death. The number of scratches, a parameter used to discriminate between browsing and grazing diets [START_REF] Semprebon | Advances in the reconstruction of ungulate eco morphology with application to early fossil equids[END_REF], is lower in the middle of the warm season than in other seasons. This implies a higher consumption of eudicots in summer, due to (1) a reduced availability of grasses in this season [START_REF] Drogoul | Nutrition et alimentation des animaux d'élevage[END_REF]; or (2) additional leguminous fodder supply. This latter strategy has been documented in Iron Age sites in nearby geographical contexts (e.g., Bouby & Rouas, 2005;[START_REF] Alonso | Using stable isotopes and functional weed ecology to explore social differences in early urban contexts: the case of Lattara in Mediterranean France[END_REF]. Another explanation (3) could be that herds were moved to nearby wooded areas, as attested by palaeoenvironmental data [START_REF] Montaner | El paleopaisatge fluvio-estuari d'Empúries[END_REF]Castanyer et al., 2016;Ejarque et al., 2022). This latter explanation could be related to a system of summer pastures.

Regarding the seasonality of animal slaughtering, cementum analysis on caprine teeth showed that the last deposit most frequently indicated the beginning of a growth zone (ca. spring) in the first period studied (550-450 BC) and the middle of a growth zone in the second phase (ca. summer). In both phases, slaughtering in winter or fall is scarce. Our results suggest that adult Caprinae were preferentially slaughtered in the spring and summer months but we must remain cautious due to the small number of samples. The lower availability of grasses and reduced grass intake in summer (Drogoul et al., 2014), confirmed by the dental microwear analysis (see section 3.5), could have prompted shepherds to slaughter adult Caprinae in spring/summer. In addition, this time frame would have been propitious to sheep shearing. For cattle, slaughtering took place all year round. However, the number of specimens analysed is even smaller than for Caprinae. These hypotheses therefore need to be tested in future studies.

CONCLUSIONS

In this paper, we presented a multiproxy approach (zooarchaeology, dental micro-and mesowear analyses, cementochronology) applied to two Iron Age sites (Empúries and Ullastret) located in the northeast Iberian Peninsula.

Dental wear analyses suggest that sheep and goats were fed in an environment dominated by shrubs and bushes in both studied periods (550-450 BC and 450-325 BC) and sites, while cattle mainly grazed in grasslands. This scenario is compatible with the hypotheses proposed in previous works: shepherds would have been forced to lead Caprinae to marginal areas less suitable for agriculture. The results are also consistent with feeding on eroded landscapes. Increased woodland clearance and land degradation in the region at that time have been documented by palaeoenvironmental studies [START_REF]Vegetation history and human activity during the last 6000 years on the central Catalan coast (northeastern Iberian Peninsula)[END_REF][START_REF]Paisatge i explotació forestal durant el I mil•leni a. n. e. a la plana empordanesa[END_REF].

In this study, cementum analyses on sheep, goat and cattle teeth provided valuable information on the seasonality of husbandry practices. Caprinae were mostly slaughtered in spring and summer. Combined cementum and microwear analyses revealed a lower consumption of grasses in summer than during the rest of the year. This could be the result of environmental/climatic factors or linked to the human control of feeding places or practices (fodder procurement).

The results obtained in this work support the hypotheses put forward in previous works. Competition between agriculture and livestock, as well as the socio-economic context, impacted livestock management and landscape use. The multiproxy approach presented here shows that animal feeding strategies, grazing areas and the season of livestock slaughtering were probably adapted to this context. This paper demonstrates how the application of this multiproxy approach can offer new perspectives and lines of research on livestock management and landscape use in Late Prehistory, which cannot be attained with traditional zooarchaeological studies.

AKNOWLEDGEMENTS

This project has been supported by LabEx AR-CHIMEDE, University Paul Valéry, Montpellier 3 ("Investissement d'Avenir" program ANR-11-

  Archaeofauna 32-1 (2023): 00-00 A total of 62 thin sections of eight sheep, 15 goat, two unidentified caprine and 13 cattle lower and upper molars were made following the protocol established by Naji et al. (2015) (Table

FIGURE 2 Mammal

 2 FIGURE 2Mammal NISP frequencies per period at Ullastret (A and B) and Empúries (C and D).

FIGURE 3 A=

 3 FIGURE 3 A= Dental mesowear results of sheep, goat and cattle; B= Microwear bivariate graph for sheep, goats and cattle from 550-450 BC. Extant data on leaf browsers (B) and grazers (G) from Solounias & Semprebon (2002). C= Microwear bivariate graph for sheep, goats and cattle from 450-325 BC. Extant data on leaf browsers (B) and grazers (G) from Solounias & Semprebon (2002). (Ull= Ullastret; Emp= Empúries; N= number; SD= standard deviation).

  NS SD SWS %XS %G N NP SD NS SD SWS %XS %G N NP SD NS SD SWS %XS

FIGURE 4 A

 4 FIGURE 4 A, B, C and D= Schematic representation of the seasonality of caprine slaughtering and teeth from Empúries and Ullastret (550-450 BC and 450-325 BC). E= Microphotograph of a longitudinal section of goat tooth showing growth zones and annuli deposits (R= resin, D= dentin, N= number).

TABLE 1

 1

	Archaeological site/Proxy	Empúries	Ullastret	Total
	Zooarchaeology	2012	1512	3524
	Dental mesowear	56	48	104
	Dental microwear	58	73	131
	Cementum analysis	19	19	38

Table 3

 3

	). Con-

TABLE 2 Mesowear

 2

score (MWS) analysis of caprines and cattle from Ullastret and Empúries (N= number of specimens; SD= standard deviation).

  Summary of dental microwear pattern of cattle, sheep and goats from Ullastret and Empúries (N= number of specimens, NP= number of pits; SD= deviation standard; NS= number of scratches; SWS= scratch width score; %XS= presence of more than four cross scratches; %G= presence of gouges).

		22 16.5 4 19.6 3.4 0.9 30	0	12 25.8 10.4 18.1 5.2 0.8 30	0	7 26.5 10.8 17.1 4	0.9 10	0
	Empúries C (5 th -4 th c. BC)	11 22.3 8.4 20.9 5.8 1.2 80	0	8 23.6 7.2 18.8 4.9 0.8 10 12.5 8 24.3 9.5 19.8 4.2 0.6 10	0
				TABLE 3

LabEx ARCHIMEDE, ANR-11-LABX-0032-01, Université Paul Valéry Montpellier 3, Montpellier, France (Received 15 February 2022; Revised 29 September 2022; Accepted 16 January 2023)ABSTRACT: Two major settlements were established in the Empordà region (north-eastern Iberian Peninsula) in the Iron Age (6 th -2 nd centuries BC), 15 km from each other: the Greek colony of Empúries and the Iberian city of Ullastret. This coexistence of two different ethnic entities in the same region -and the same environment -presents an opportunity to shed light on the variability of livestock practices in these settlements, as well as their relationship with the environment. For this purpose, in addition to traditional zooarchaeological approaches, we reconstructed the animal palaeodiet and seasonality of death using three different proxies: dental meso-and microwear analyses, and cementum analysis. The results of the study, conducted on caprine and bovine teeth from both sites, support the hypothesis proposed in previous works, namely that the environment was exploited differently depending on livestock species. Caprinae (sheep and goats) were fed in marginal areas and less suitable areas for agriculture while, in contrast, cattle grazed on rich grassland and/or in wetland areas. The combination of cementochronology and dental microwear suggests a lower consumption of grasses in summer for Caprinae and cattle. Seasonality results show that Caprinae were mainly slaughtered during spring and summer.

Archaeofauna 32-1 (2023): 00-00 LABX-0032-01), and a Grant FJC 2021-046527-I funded by MCIN/AEI/10.13039/501100011033 and "European Union NextGenerationEU/PRTR".

Archaeofauna 32-1 (2023): 00-00 SUPPLEMENTARY MATERIAL