Clément Perrin
email: perrin@ipgp.fr

Felix Waldhauser

Eunseo Choi

Christopher H Scholz

California, story. Clément Parkfield

Perrin

Persistent fine-scale fault structure and rupture development: A new twist in the Parkfield, California, story

Introduction

The Parkfield section of the San Andreas fault is one of the best-instrumented earthquake source region in the world. Over the past 50 years, one and a half seismic cycles, including two M~6 earthquakes in 1966 and 2004, have been observed there with full geodetic and seismic coverage. Yet it is still unclear how the earthquake generation process takes place in this peculiar area.It has long been recognized that the M~6 earthquakes break a stuck patch imbedded in an otherwise creeping fault -first noted by [START_REF] Harris | Detection of a Locked Zone at Depth on the Parkfield, California, Segment of the San Andreas Fault[END_REF] and most recently modeled by [START_REF] Barbot | Under the Hood of the Earthquake Machine : Toward Predictive Modeling[END_REF]. It has also been observed that the rupture zone lies between two geometrical discontinuities in the fault, a slight compressional bend in the north near Middle Mountain (MM) and an extensional jog in the south near Gold Hill (GH) [START_REF] Lindh | Control of Rupture by Fault Geometry During the 1966 Parkfield Earthquake[END_REF]. Neither the connection between these features nor their role in the nucleation and termination of the M~6 earthquake ruptures have been elucidated.

The 1966 earthquake initiated at the northern end of the rupture area near Middle Mountain, propagated unilaterally to the southeast along the entire Parkfield section and stopped near Cholame (CH) (Fig. 1 and 2; [START_REF] Brown | The Parkfield-Cholame, California, earthquakes of June-August 1966 -surface geologic effects, water-resources aspects, and preliminary seismic data[END_REF][START_REF] Mcevilly | The Parkfield, California, earthquakes of[END_REF].

In contrast, the 2004 earthquake initiated in the south near Gold-Hill (GH), then propagated mainly to the northwest and stopped near MM (Fig. 1 and 2; [START_REF] Langbein | Preliminary Report on the 28 September 2004, M 6.0 Parkfield, California Earthquake[END_REF][START_REF] Rymer | Surface fault slip associated with the 2004 Parkfield, California, earthquake[END_REF]. Both surface ruptures were similar, formed by small fractures associated with local displacements less than 0.1 m [START_REF] Brown | The Parkfield-Cholame, California, earthquakes of June-August 1966 -surface geologic effects, water-resources aspects, and preliminary seismic data[END_REF][START_REF] Rymer | Surface fault slip associated with the 2004 Parkfield, California, earthquake[END_REF]. Source inversion models of the 2004 event reveal two main slip patches close to either end of the rupture extent, with maximum slip ranging from 0.5 to 1m (see an example in Fig. 3), followed by post-seismic afterslip near the surface above the coseismic rupture zone (e.g., [START_REF] Barbot | Postseismic deformation due to the Mw 6.0 2004 Parkfield earthquake: Stress-driven creep on a fault with spatially variable rateand-state friction parameters[END_REF][START_REF] Johanson | Coseismic and postseismic slip of the 2004 Parkfield earthquake from space-geodetic data[END_REF][START_REF] Langbein | Coseismic and initial postseismic deformation from the 2004 Parkfield, California, earthquake, observed by global positioning system, electronic distance meter, creepmeters, and borehole strainmeters[END_REF][START_REF] Murray | Slip on the San Andreas fault at Parkfield, California, over two earthquake cycles, and the implications for seismic hazard[END_REF]) but also at depth in the lower crust [START_REF] Bruhat | Evidence for postseismic deformation of the lower crust following the 2004 Mw6.0 Parkfield earthquake[END_REF]. Similar co-seismic distribution associated with shallow postseismic behavior has been observed for the 1966 earthquake (see an example in Fig. 3; e.g. [START_REF] Archuleta | Dynamic rupture in a layered medium: the 1966 Parkfield earthquake[END_REF][START_REF] Scholz | Seismic and aseismic slip on the San Andreas Fault[END_REF][START_REF] Segall | How similar were the 1934 and 1966 Parkfield earthquakes?[END_REF].

In order to better understand rupture initiation, propagation and arrest of the two latest M6 earthquakes at Parkfield, we present here a combined analysis of high-resolution near surface and earthquake data that image the fine-scale structures of the San Andreas Fault. We bring together geologic fault information available in the literature and a broad range of satellite images and topographic data to refine the long-term surface fault traces along the Parkfield section. We also compute high-resolution seismicity recorded at the Northern California Seismic Network (NCSN) between 1969 and 2015, and aftershocks recorded during the 2 months following the 1966 earthquake [START_REF] Eaton | Aftershocks of the 1966 Parkfield-Cholame, California, earthquake: A detailed study[END_REF], spanning a nearly 50 year time period including two ~M6 events. Using these new data, we relate the structures observed at the surface to persistent deformation patterns along the three-dimensional fault plane at depth. Supported by numerical simulations of the long-term fault deformation, we propose a new mechanical model for the Parkfield segment of the San Andreas fault that differs from previous studies (e.g., [START_REF] Simpson | San Andreas fault geometry in the Parkfield, California, region[END_REF].

Data analysis and results

Surface fault traces

According to the quaternary surface fault map from the USGS (U.S. Geological Survey and California Geological Survey, 2006); blue lines Fig. 1b and1c), the ~30 km long Parkfield section is composed of two linear, roughly N139°E striking, sub-parallel strands (also called the main San Andreas fault and the southwest fracture zone) that join each other at both terminations of the Parkfield section through two structural complexities : i) at the southern end near GH, an extensional jog, where the fault trace bends twice with deviations up to 28° from the mean fault strike (GH-CH, Fig. 1b and2c); ii) at the northern end, near MM, a 5° restraining bend [START_REF] Lindh | Control of Rupture by Fault Geometry During the 1966 Parkfield Earthquake[END_REF] (MM, Fig. 1b and2c).

We analyzed multi-resolution satellite images (Google Earth, Landsat 7) and topographic data (SRTM 30 m, LiDAR available from opentopography.org) and identified additional structural features (black lines in Fig. 1b and1c): around MM, the active fault trace is surrounded by multiple secondary faults several kilometers away from the main San Andreas Fault and clearly imprinted in the topography (Fig. 1c). Most of them have been previously recognized as secondary strike-slip and reverse faults in the field [START_REF] Thayer | Structural Geology of the San Andreas Fault Zone at Middle Mountain, near Parkfield, Central California[END_REF]. Altogether they contribute to forming the localized relief observed at the center of the main fault trace. The overall shape of this structural feature, which encompasses an area approximately 5 km wide and 10 km long, resembles a restraining step-over feature that separates the locked Parkfield section from the creeping section.

These secondary faults strike between N110°E and N150°E, in good agreement with earthquake locations derived in this study and focal mechanisms from Hardebeck (2010) that show oblique reverse and normal faulting (Fig. 1c; see also [START_REF] Thurber | Three-Dimensional Compressional Wavespeed Model, Earthquake Relocations, and Focal Mechanisms for the Parkfield, California, Region[END_REF]. The 1966 and 2004 earthquakes initiated and stopped near these large structural discontinuities observed at the surface near GH-CH and MM.

Earthquake relocation

We use seismicity data to investigate the depth extent of the heterogeneities observed near the surface. We first apply a 3D grid search method (NLLoc program; [START_REF] Lomax | Probabilistic earthquake location in 3D and layered models: Introduction of a Metropolis-Gibbs method and comparison with linear locations[END_REF] to determine new absolute locations of more than ~17,000 earthquakes from 1966-2015 in the 3D velocity model of [START_REF] Thurber | Three-Dimensional Compressional Wavespeed Model, Earthquake Relocations, and Focal Mechanisms for the Parkfield, California, Region[END_REF]. The strike of the new 3D absolute locations is rotated slightly (up to ~2°) counterclockwise compared to the original NCSN locations as a result of properly accounting for the velocity variations across the fault zone. We then use the double-difference algorithm HypoDD [START_REF] Waldhauser | HypoDD: A computer program to compute double-difference hypocenter locations[END_REF][START_REF] Waldhauser | A Double-difference Earthquake location algorithm: Method and application to the Northern Hayward Fault, California[END_REF] to compute precise relative locations using both travel-time differences from P-wave pick data (available from 1966 to 2015) and P-and S-wave cross correlation delay time data [START_REF] Schaff | Waveform cross-correlation-based differential traveltime measurements at the northern California seismic network[END_REF]. We included all events with at least 4 stations in order to be able to simultaneously relocate the aftershocks of the 1966 earthquake [START_REF] Eaton | Aftershocks of the 1966 Parkfield-Cholame, California, earthquake: A detailed study[END_REF] together with the rest of the seismicity. The resulting locations are shown in Fig 1,S1c and S1d. These locations resolve the fine shape of the fault plane at a level of a few tens of meters. The median lateral relative error at the 95% confidence level, derived from a bootstrap analysis of the fine residuals, is 0.039 km and the median vertical relative error is 0.022 km (Fig. S2).

Most events with M < 1 (ML is considered for events M > 3, MD for smaller events; see [START_REF] Uhrhammer | California integrated seismic network (CISN) local magnitude determination in California and vicinity[END_REF] are located in the creeping section (see Supplementary Fig. S1a and S1b). A cross-section at MM (Fig. 1d) indicates that some events also delineate secondary oblique strands in the eastern part (see black arrows) that connect to the main fault at ~10 and ~12 km depth, close to the 1966 nucleation area. Several reverse and normal faulting focal mechanisms located between 10 and 15 km deep (events 3, 4, 5 and 6 in Fig. 1c and1d) are aligned with those secondary strands, while those located between 0 and 10 km deep (events 1, 2 and 7 in Fig. 1c and1d), appear to be related to the secondary fault traces observed at the surface (black crosses, Fig. 1d).

Distribution of the 1966 and 2004 aftershocks

In the following analysis, we only consider earthquakes of M ≥ 1, corresponding to the minimum magnitude of completeness of the earthquake catalog (see Supplementary Fig. S3). Fig. 2 shows the distribution of the 1966 and 2004 aftershocks within 2 months of the mainshock (in blue and red, respectively). The along strike distribution of the aftershocks is sharp and fairly similar for both events, mostly situated between 4 and 12 km deep along the Parkfield section. Distinctive seismicity patches (e.g. streaks, clusters, repeating events, [START_REF] Waldhauser | Streaks, multiplets, and holes: High-resolution spatio-temporal behavior of Parkfield seismicity[END_REF][START_REF] Thurber | Three-Dimensional Compressional Wavespeed Model, Earthquake Relocations, and Focal Mechanisms for the Parkfield, California, Region[END_REF] outline an area devoid of earthquakes between 6-9 km deep, corresponding to the main locked patch of the fault that broke in the 1966 and 2004 events (dashed lines in Fig. 2b, [START_REF] Harris | Detection of a Locked Zone at Depth on the Parkfield, California, Segment of the San Andreas Fault[END_REF].

The 1966 and 2004 aftershock sequences are similarly distributed, especially near GH-CH where two distinct patches of seismicity, at ~5 km and ~10 km depth, are superimposed (Fig. 2b), suggesting that the same main fault strand has been activated during the two earthquakes. A high concentration of number of aftershocks (including all M ≥ 4 events) and cumulative moment release are observed near the rupture tips of the 1966 and 2004 earthquakes (Fig. 2a,b; Fig. 3a,b).

We note that in the northwestern part of the Parkfield section, the 1966 aftershocks tend to deviate further to the northeast away from the San Andreas fault trace (blue dots in map view on Fig. 2a). This deviation is also observed in the original absolute locations and might be due to the poor station coverage in the northwestern part of the 1966 rupture zone [START_REF] Eaton | Aftershocks of the 1966 Parkfield-Cholame, California, earthquake: A detailed study[END_REF]. However, few events are concerned and our relocations are in good agreement with previous studies in the southeastern part of the Parkfield section, where most of the seismicity is localized (Fig. 2; e.g., [START_REF] Thurber | Three-Dimensional Compressional Wavespeed Model, Earthquake Relocations, and Focal Mechanisms for the Parkfield, California, Region[END_REF][START_REF] Waldhauser | Large-scale relocation of two decades of Northern[END_REF]. More generally, less than 20 permanent and temporary stations recorded aftershock activity following the 1966 earthquake, with several periods of interrupted service [START_REF] Eaton | Aftershocks of the 1966 Parkfield-Cholame, California, earthquake: A detailed study[END_REF]. This likely explains why the aftershock activity and the cumulative seismic moment release shown in Figure 3 are much lower for the 1966 sequence compared to the one in 2004. 1969-2004 and 2005-2015 We choose to represent in Fig. 2 only the best (i.e., correlation based) relocations (gray dots) for the two interseismic periods from 1984-2004 and 2005-2015 for a fine-scale representation of the geometrical features (see Supplementary Fig. S1 for the entire relocated catalog). The interseismic seismicity distribution is remarkably similar to the co-and postseismic distributions. Figure 3 shows along strike cross sections of the earthquake density and cumulative seismic moment released within 2 months after the 1966 and 2004 Parkfield earthquakes (Fig. 3a and3b, respectively), and during the [START_REF] Stein | Seismic potential revealed by surface folding: 1983 Coalinga, California, earthquake[END_REF][START_REF] Waldhauser | Streaks, multiplets, and holes: High-resolution spatio-temporal behavior of Parkfield seismicity[END_REF]and 2005-2015 interseismic periods (Fig. 3c and 3d, respectively) interseismic periods (Fig. 3c and 3d, respectively). Density contour lines point out five, possibly six, persistent kilometer scale seismicity patches during the interseismic periods (Fig. 3c and3d). A patch of higher density is located at the northern end of the Parkfield section at 4 km depth, a patch that was also activated by aftershocks of the 2004 event (Fig. 3b).

Interseismic periods

The interseismic, cumulative moment release maps allow us to distinguish other patches which are lining up at 4-6 km and ~10 km depth (Fig. 3c and3d), again in good agreement with patches that were active during the 1966 and 2004 post seismic periods (Fig. 3a and3b).

3D fault geometry

We determine the three-dimensional fault geometry by applying a principal component analysis to all events within 3 x 3 km wide boxes (in map view), stepping at 1 km intervals along the fault trace. For each box, we calculate strike and dip of the plane that minimizes the distance between hypocenters and fault surface during each of two interseismic periods (1984-2004 and 2005-2015) in an attempt to capture the longterm signature of the fault geometry (Fig. 2c). For simplicity, we assume, in each box, a constant dip of the calculated plane. The calculated strike and dip values are presented in Fig. 2c, while an animated view of the data and the combination of the best fitting planes is shown in Supplementary Fig. S6.

The strike measurements show that the strike along the mainshock rupture zone, which includes the stuck patch, is fairly constant, while it deviates from that trend at both ends: the strongest deviations range from N125°E to N150°E (at MM and GH), consistent with surface measurements (see black line in Fig. 2c), albeit with fluctuating amplitudes. The overall strike variations are also in good agreement with focal mechanisms derived by [START_REF] Hardebeck | Seismotectonics and fault structure of the California Central Coast[END_REF]Supplementary Fig. S7). Figure 2d shows the standard deviations of the seismicity from the best fitting planes for each interseismic periods, characterizing the normal distance of the events from the best fitting planes (see also Supplementary Furthermore, the fault dip varies steadily from northeast dipping in the north to southwest dipping in the south, indicating a twist in the fault plane along a horizontal axis (Fig. 2c; e.g. [START_REF] Kim | Seismicity and fault geometry of the San Andreas fault around Parkfield, California and their implications[END_REF]. The 3D fault geometry of the San Andreas fault at Parkfield can thus be best described as a helicoid (Supplementary Fig. S6). The change in fault dip from northwest to southeast can also be observed in the original NCSN locations, as well as the 3D grid-search absolute locations, indicating that the twist is not an artifact of the relocation procedure or the model used to locate the events (see Supplementary Fig. S8).

Because the seismicity patches are well defined and aligned, we used the principal component analysis to minimize the orthogonal hypocenter distances to a single fault plane in each box (Supplementary Fig. S4). However, we also tried to investigate the possibility of multiple fault strands or variations of the fault dip at depth below and above 7 km (i.e. the limit between the main patches of seismicity; see Supplementary Fig. S4). The resulting fault planes were poorly constrained, especially the dip values, because their fits were mostly controlled by single sub-horizontal streaks, rather than events well distributed on a plane in three-dimensional space. Estimating an orientation by fitting one single plane in each box tend to smooth local irregularities that might be due to multiple fault strands at depth. Event distributions away from the fault planes (Supplementary Fig. S5) show that the Parkfield section is mainly characterized by a clear and sharp single strand of ~100m wide. Secondary peaks, indicating possible oblique secondary strands, are observed around MM and GH-CH. Therefore, our approach cannot describe in detail the complexity of the damage zone around the fault core, but it allows us to highlight significant discontinuities with great confidence.

Numerical model of the effect of a strong stuck patch on fault geometry

The creeping section is known from studies at SAFOD to be very weak, with a friction coefficient µ~0.15 due to the presence of weak phyllosilicates in the fault zone (i.e. saponite at shallower depth, [START_REF] Carpenter | Frictional properties of the active San Andreas Fault at SAFOD: Implications for fault strength and slip behavior[END_REF][START_REF] Lockner | Low strength of deep San Andreas fault gouge from SAFOD core[END_REF], to talc below 3 km deep, Moore et al., 2016). They appear to be produced by reaction between the serpentinite of the Coast Ranges ophiolite and the silicious rocks of the Salinian block [START_REF] Moore | Chemical controls on fault behavior: Weakening of serpentinite sheared against quartz-bearing rocks and its significance for fault creep in the San Andreas system[END_REF]) -a geologic context that is only present along the creeping section [START_REF] Allen | The tectonic environments of seismically active and inactive areas along the San Andreas fault system[END_REF][START_REF] Irwin | Effects of geological structure and metamorphic fluids on seismic behavior of the San Andreas fault system in central and northern California[END_REF]. On the other hand, the stuck patch at Parkfield cannot be lined with weak phyllosilicates because they are velocity strengthening and will not facilitate seismic instability. The only possible materials that can be velocity weakening and exhibit the frictional strength of the stuck patch in this area are granite, antigorite or lizardite forms of serpentinite [START_REF] Moore | Talc-bearing serpentinite and the creeping section of the San Andreas fault[END_REF]. These all have high strength with friction coefficients between 0.5-0.75. Thus, the seismogenic stuck patch has a higher friction coefficient than the surrounding creeping segment and thus can sustain greater shear stress than the stress-drops observed in Parkfield earthquakes (on the order of ~1 MPa). Since the earthquakes do not relieve all the stresses on the stuck patch, there must be a permanent torque produced by that high-stress patch surrounded by low stress regions. As a result, over time periods longer than the seismic cycle, slow inelastic deformation will occur in response to this torque and the twist in the fault plane will develop.

In the following section, we present a numerical model that supports the above idea.

Rather than simulate the entire process of accumulating torque and inelastic deformations, we focus on showing that a high-friction patch embedded in a low-friction fault plane, representing the stuck patch situation, is able to generate deformations that are consistent with the proposed helicoidal twisting of the fault plane.

Model setup

We use PyLith, a portable, scalable finite element software for simulation of crustal deformation and earthquake faulting [START_REF] Aagaard | PyLith User Manual[END_REF][START_REF] Aagaard | A domain decomposition approach to implementing fault slip in finite-element models of quasi-static and dynamic crustal deformation[END_REF], to simulate the deformation of linear elastic crust containing a straight strike-slip fault with heterogeneous frictional strengths. The model domain, covering an area 40 × 162 km wide and 21 km-thick (Fig. 4a), is discretized into hexahedral linear elements with edge length of 1 to 1.5 km. Density of the crust is assumed to be 2500 kg/m 3 and Vs and Vp to be 3 and 5.2915 km/s. These values are equivalent to shear and bulk modulus of 22.5 and 55 GPa, respectively. Completely contained within the crust, the modeled zerothickness fault plane is 90 km long and 20 km tall, going through the center of the domain parallel to the longest edge of the domain. The creeping section is represented by a low-friction section with a uniform friction coefficient of 0.1 and a cohesion of 2 MPa ("creeping fault plane" in Fig. 4a). The stuck patch is modeled as a simplified rectangular zone of 20 × 5 km ("high friction patch" in Fig. 4; [START_REF] Harris | Detection of a Locked Zone at Depth on the Parkfield, California, Segment of the San Andreas Fault[END_REF]. It has a uniform friction coefficient of 0.6 and a cohesion of 2 MPa. The fault section that broke during the 1857 earthquake ("high-friction fault plane" in Fig. 4), is considered to have the same friction coefficient and cohesion as the stuck patch. All the fault sections are allowed to slip because we model fault behaviors for a time period corresponding to multiple earthquake cycles, over which no part of the Parkfield section would remain completely locked.

Right-lateral motion at a rate of 2 cm/yr on two fault-parallel domain boundaries (relative on-fault slip rate of 4 cm/yr) with the normal component of velocity (vx) set to be zero and the vertical component (vz) free (Fig. 4b). On the fault-perpendicular sides, a simple-shear velocity field is applied such that the strike-parallel component varies between +2 and -2 cm/yr linearly with distance across the fault. The fault perpendicular component is set to be zero and the vertical component is free on these boundaries. The bottom boundary is a free-slip surface. The top boundary is a free surface with zero traction. The effect of gravity is considered and the initial stress is set to be lithostatic.

We get fault displacements and stresses as a response to incrementally increasing amounts of boundary displacements. Although we apply velocities to generate strikeslip kinematics, the velocity boundary conditions are simply a way of increasing strikeslip motions incrementally at each time step. Our models simulate only linear elastic behaviors of rocks and does not directly show rate-dependent effects. Nonetheless, the correspondence principle (e.g., Sec. 2.5.7, LeMaitre, 2001) ensure the relevance to permanent helicoidal twisting of the Parkfield section. According to this principle, strain rate field in a linear viscoelastic model would be equivalent to elastic strain field in the corresponding linear elastic model under the same boundary conditions. So, we can view our model's strain field as a strain rate field for long term deformation.

Results

Figure 5 presents the northeastern displacements on the fault plane when the highfriction patch and the upper 2 km of the high-friction plane started slipping after 112 m of cumulative displacement (Fig. 5a and5c). After an unspecifiable number of earthquake cycles, slip would accumulate on the high-friction patch as well. Thus, we also look at the northeastern displacements when the entire high-friction patch as well as the upper 5 km of the high-friction plane has slipped after 160 m of boundary displacement (Fig. 5b and5d). Relative to the state in Fig. 5a and5c, northeastern displacements have increased at the northern tip and more at the southern tip of the patch. A major warp of the fault plane occurs at the northern tip of the high-friction patch, corresponding to MM, with northeastern displacements > 1 m along a zone approximately 10 km long (Fig. 5a and5c). We note also that a section of the fault situated in the southern tip of the stuck patch is slightly displaced toward the southwest when most of the high-friction patch is locked (Fig. 5a and5c), accentuating the general warp of the fault plane. When the high friction patch starts to slip, the whole section is progressively shifted toward the northeast, but the general warping of the fault plane is still preserved (figure 5b and5d). An animated view of the model is available in the Supplementary material (see Supplementary Fig. S9).

Cross sections in Fig. 6 show the deformation of the San Andreas Fault plane at specific locations along the model (Fig. 6a) after the stuck patch has been fully mobilized (i.e., the stage in Fig. 5b and5d). The vertical cross-section at MM (Fig. 6b) indicates a northeastern displacement of the fault plane from its original location near the surface (hence a southwest dipping fault). Most of the deformation diffused in the crustal medium is fairly symmetrical relative to the fault plane. It seems to root at the base of the stuck patch and to widen towards the surface in a flower-like structure. The deformation is less pronounced at the southern tip of the stuck patch (i.e., GH, Fig. 6c).

Yet, we can distinguish a slight deflection of the fault plane, which is displaced toward the northeast in the shallower and deeper part of the model (i.e. the creeping zones) and less so in an area that corresponds to the stuck patch. Thus, the dip of the fault plane changes with depth (i.e. toward the southwest, northeast and southwest). Figure 6d is a horizontal cross-section situated at 7.5 km depth, crossing through the creeping zones and the stuck patch (see Fig. 6a). It highlights the variation of the fault strike (grey dotted line) compared to a referenced linear fault direction (black dotted line), in the same way as in Fig. 2d. As shown in Fig. 5, the fault plane is mostly deformed toward the northeast at MM and GH, and much less so in the southern half part of the stuck patch.

The strike variations within the stuck patch are consistent with the clockwise rotation of the stuck patch inferred from the seismic data (Fig. 2c,2d).

Figure 7 provides additional information on on-fault stress distribution at two stages: after 112m of displacement (Fig. 7a), the maximum shear stress (~100 MPa) occurs at the left (northwestern) tip and corners of the rectangular stuck patch (corresponding to MM), which begins to slip. At this point, slip along the stuck patch, which is meant to represent cumulative slip from repeated earthquakes, is about 30 cm with a stress drop of ~ 1 MPa, which is regained in ~ 30 yrs. Yet, after 160 m of cumulative displacement, stress remains at about 100 MPa throughout the seismic cycles, even after the highfriction patch has slipped (Fig. 7b). Thus, slow anelastic processes can produce permanent deformation of the same type as the elastic ones shown in Fig 5 . This result indicates that the high-friction patch can uphold a higher level of stress than the surrounding creeping section permanently, i.e., over a long time period corresponding to many earthquake cycles. In both states, accumulating slip on the creeping section generates high shear stress near the clamped bottom edge of the fault plane. Likewise, shear stress is concentrated right beneath the slipping region on the high-friction section, southeast of the patch (Fig. 7b). Profiles of shear and normal stresses (σxy and σxx, respectively) are shown in Supplementary Fig. S10. Along dip, the shear stress increases on the stuck patch relative to the surrounding creeping regions and higher stress concentrations are localized at the lower part of the stuck patch. Along strike, the variation of shear and normal stresses are consistent with compression at the left tip (MM) and extension at the right tip (GH-CH) of the stuck patch, resulting in clockwise rotation.

Discussion

Progressive deformation of the San Andreas fault over the long-term.

Shear stress and friction contrast at the origin of the twist

Our observational results show that the San Andreas fault plane within the rupture region of the recent ~M6 events has been deformed into a helicoidal surface that extends from the compressional bend at MM in the north to the extensional jog near GH in the south (Fig. 2c and Supplementary Fig. S6). It connects the northern part of the fault to the Cholame section in the south through a ~ 2 km wide step. By incorporating changes in fault-slip behavior at depth, our numerical model (Fig. 6) reproduces this twist and thus the variation of the fault strike (Fig. 2d,6d) and fault dip (Fig. 2d,6b,6c) deduced from the seismicity analysis. The stuck patch is stronger than the surrounding creeping areas by a factor of ~4 (see section 3). The contrasting high right-lateral shear stresses on the patch will induce a torque that will bend the stuck section of the fault clockwise, as shown in Fig. 8a (map view). Because the shear stresses increase with depth, another torque increases the clockwise bending with depth (Fig. 8b). In our numerical model, the increase in shear stress from the surface to 10 km depth is about 30 MPa in the creeping section away from the stuck patch (Fig. 7a,b). However, a much stronger torque is resulted, about 120 MPa over the same vertical distance, through the stuck patch (Fig. 7b). This result shows how the assumed contrast in frictional strength can be responsible for the torques localized around the stuck patch. The two mechanisms in Fig. 8a and Fig. 8b would produce the helicoidal fault geometry that we observe (Fig. 8c). This deformed section acts like a hinge and gate, with the hinge in the north marked by the compressional bend at its onset (MM in pink, fig. 8c), and the open end of the gate to the south, where it is accommodated by the extensional jog (GH-CH in blue, fig. 8c). This suggests that the geometrical complexities at the ends of the stuck patch result from the contrast of strength of the stuck patch with the surrounding creeping section of the fault.

Formation and evolution of geometrical complexities and secondary structures.

One hypothesis would consider that the geometrical complexities are inherited structures formed before the presence of the stuck patch in the geologic evolution of the San Andreas fault. They would correspond to ancient relay zones between fault segments, where many WNW-ESE long-term secondary structures (Fig. 1c) would be connected to the main San Andreas Fault. As such they could play an important role in the fault structure and the rupture development of ~M6 Parkfield earthquakes. They would contribute to control the extent of the rupture, acting as structural barriers (e.g., [START_REF] Biasi | Bends and Ends of Surface Ruptures[END_REF][START_REF] Manighetti | Earthquake scaling, fault segmentation, and structural maturity[END_REF][START_REF] Stirling | Fault trace complexity, cumulative slip, and the shape of the magnitude-frequency distribution for strike-slip faults: a global survey[END_REF] and places where stress concentrations might also favor rupture initiation (e.g., [START_REF] Aki | Characterization of barriers on an earthquake fault[END_REF][START_REF] Barka | Strike-slip fault geometry in Turkey[END_REF][START_REF] Manighetti | Generic along-strike segmentation of Afar normal faults, East Africa: Implications on fault growth and stress heterogeneity on seismogenic fault planes[END_REF][START_REF] Shaw | Initiation propagation and termination of elastodynamic ruptures associated with segmentation of faults and shaking hazard[END_REF]. The secondary structures are possibly oblique thrusts and folds formed in a transpressive regime during the southern propagation of the San Andreas Fault (e.g. [START_REF] Perrin | Off-fault tip splay networks: A genetic and generic property of faults indicative of their long-term propagation[END_REF]Titus et al., 2011 and references therein). Their obliquity varies depending on the distance from the main San Andreas Fault [START_REF] Titus | Geologic versus geodetic deformation adjacent to the San Andreas fault, central California[END_REF] and some of them could be still potentially active (ex: the 1983, Mw 6.2 Coalinga thrust earthquake, [START_REF] Stein | Seismic potential revealed by surface folding: 1983 Coalinga, California, earthquake[END_REF][START_REF] Das | An elliptical asperity in shear: fracture process and seismic radiation[END_REF], Mw 6.1 Kettleman Hills thrust earthquake, Ekström et al., 1992). The diffuse seismicity around Middle Mountain (Fig. 1, 2, S1, S4 and S5) might be explained by the presence of these secondary faults in the crustal medium surrounding the main San Andreas Fault. Yet, it seems unlikely that the presence of the stuck patch randomly coincides or results from the presence of the structural complexities. The San Andreas fault is a large cumulative slip fault with a very smooth (i.e., linear) surface fault trace along the central section (Wesnousky, 1988), except at the tips of the Parkfield section. Thus, it is necessary to find another mechanism responsible for the local twist of the fault plane. We propose a second, more likely, hypothesis: the local structural complexities were lately formed by the torque of the fault plane due to the presence of a locked fault patch surrounded by creeping zones at depth. Thus, secondary geological fault traces at MM and GH-CH should be considered in the common view of the active fault traces in California. They reflect the stress concentrations at the tip of the stuck patch at depth but can also enhance them. The geometrical complexities are marked by a wider distribution of the seismicity (Fig. 2d and Supplementary Fig. S4 andS5), associated with a deflection of the San Andreas fault plane (Fig. 8 and Supplementary Fig. S6). Finally, the flower structure shape highlighted in our model (Fig. 6b) is in good agreement with the development of secondary structures at depth that widen towards the surface (Fig. 1d and Supplementary Fig. S6). [START_REF] Simpson | San Andreas fault geometry in the Parkfield, California, region[END_REF], using geological and geodetical observations in addition to the 2004 aftershock locations of [START_REF] Thurber | Three-Dimensional Compressional Wavespeed Model, Earthquake Relocations, and Focal Mechanisms for the Parkfield, California, Region[END_REF], interpreted a deflection of the active San Andreas fault near GH-CH along the vertical axis in the upper ~6 km from its surface expression as the result of non-elastic behavior of upper crustal rock units. Our observation of a systematic change in fault dip along the entire Parkfield section from MM to GH-CH suggests an additional torque along the horizontal axis that causes the helicoidal nature of the three-dimensional fault plane, as shown in our model (Fig. 6). This shape is not predicted by the 2D numerical model in [START_REF] Simpson | San Andreas fault geometry in the Parkfield, California, region[END_REF], which considers a constant, near-vertical dip at depth to obtain a warp in the fault plane with a major deflection at GH-CH, but not at MM. They suggest that the warp is due to slip in the creeping zone (including the Parkfield stuck patch), which results in inelastic deformation that is not released further to the south by large events similar to the 1857 earthquake. Our model of the helicoidal twist is based on the long-term difference in strength between persistent structures that represent stuck and creeping areas, and hence does not require inelasticity with such a short time constant.

Interpretations and numerical modeling of the twist at Parkfield: strengths and limitations

The high friction patch (µ ~0.6) is velocity weakening and so exhibits an earthquake cycle. It is embedded in a region of low friction (µ ~0.2) which continuously slides in a stable manner. This strength contrast will place a clockwise torque on the system which by slow inelastic deformation over geologic time has twisted the fault. In this regard, our numerical models serve only two purposes: i) to demonstrate that the clockwise torque can arise from the assumed contrast in frictional strength along the Parkfield section and ii) to show that the clockwise torque can be sustained over many earthquake cycles on the high-friction patch.

In our model, the shear stress magnitude is significant (~100 MPa) and the correspondend principle between the linear elastic and the linear viscoelastic strain implies, over the long term, that our model's strain field could induce the development of permanent inelastic deformation in the upper part of the crust. This deformation might correspond to the secondary faults and fractures observed at surface and at depth at the tips of the stuck patch. More sophisticated fault and material models would be needed to investigate a causal relationship between the helicoidal twist of the stuck patch and the formation of secondary structures.

Our model indicates lateral deformation of the San Andreas fault of 1.5 m to 2 m at MM and 1.3 m to 2.5 m at GH (at surface and at depth, respectively) after 160 m of cumulative slip. Considering the current total cumulative slip around Parkfield since the fault initiation (i.e. ~315 km; [START_REF] Crowell | The San Andreas fault system through time[END_REF]Revenaugh and Reasoner, 1997 and references therein), it would imply a total northeastern fault deformation of ~3 to ~4 km and ~2.5 to ~5 km at MM and GH, respectively. While this deformation is comparable observations at GH (~2 km wide step-over at the surface), it is not the case at MM. However, we note that our model has exaggerated contrast in friction strength between the patch and the creeping section. Increasing the friction coefficient of the creeping section to 0.2, compared to 0.1 used here, would possibly subdue the amount of the northeastern fault deformation. In addition, our model does not include the development of new secondary faults and fractures, which may have been created in the surrounding medium to accommodate the deflection of the San Andreas fault, just like the ~5 km wide zone of deformation identified at MM (see Fig. 1 and Thayer, 2006).

Additional minor faults might be also present but they have been possibly eroded or hidden under sedimentary deposits.

Finally, we remind that our model is a simplified view that cannot exactly reproduce the deformation pattern of the San Andreas fault at Parkfield. The results are sensitive to the initial conditions (i.e. relative location of the stuck patch at depth and connection with the locked section in the south). For instance, we assume a vertical initial fault dip and the model does not reproduce the NE dip of the fault in the creeping section (Fig. 2c). On the other hand, we are not sure that this dip is related to the torque of the fault plane due to the presence of the stuck patch or rather to the initial fault dip of the San Andreas fault. All these parameters were not explored in this study and need to be tested in future models. Our model is used here to show that the right lateral motion and the change in fault sip behavior due to the presence of the stuck patch are in good agreement with the torque highlighted in our observations. The resulting stress and displacement values should be used with caution.

Impact of fault structures on the seismicity behavior and the earthquake cycle

Our results confirm that major structural discontinuities, such as secondary reverse and normal faults, are present at the surface near the transitions between the Parkfield and creeping sections at MM, and the Parkfield and Cholame sections near GH-CH (Fig. 1).

Fine-scale analysis of the seismicity indicates that these structural complexities extend at depth and correlate with the vertices of the elliptical stuck patch that ruptures in M6 events (Fig. 2 and3). It is recognized that aftershocks occur mostly at the edges of regions experiencing high coseismic slip or in areas of low coseismic slip (e.g., [START_REF] Mendoza | Aftershock patterns and main shock faulting[END_REF], regardless of initial stress conditions. Along the Parkfield section, earthquake locations in both the post-seismic and interseismic periods over one and a half seismic cycles delineate the margins of the stuck patch (i.e. transition between locked and creeping areas on the fault; Fig. 2 and3). Especially, figure 3d updates the common view of the "stationary" seismicity around the stuck patch, going through a new interseismic period, during the decade following the 2004 earthquake. Therefore, our observations indicate that the stress concentrations that surround the stuck patch and therefore the stuck patch itself are long-term permanent features controlling the rupture process.

Based on the results from our model (Fig. 7), we schematize in Fig. 9 the state of stress and seismicity distribution on the fault plane. The seismicity distribution reflects significant stress concentrations between creeping regions, which are velocity strengthening and hence frictionally stable, and the locked region (i.e. stuck patch), which is velocity weakening and hence unstable. During interseismic periods (stage 1 in Fig. 9), the stress increases on the stuck patch relative to the surrounding creeping regions and higher stress concentrations are localized at the edge of the stuck patch. The resulting stress concentrations lead to higher seismic activity at the transition between locked and creeping regions. Coseismic slip produces surrounding stress concentrations during the postseismic period (stage 2 in Fig. 9) which also induces enhanced seismicity in that area. Then, the stress accumulates again along the locked section during the following interseismic period (stage 3 in Fig. 9), producing a similar behavior as in stage 1.

Our model shows that the long-term deformation of the fault plane is prominent at the along-strike extremities of the stuck patch (Fig. 5 and6), leading to a concentration of stress and strain that can affect the seismic cycle. These zones also correspond to the nucleation locations of the 1934[START_REF] Mcevilly | The Parkfield, California, earthquakes of[END_REF][START_REF] Waldhauser | Streaks, multiplets, and holes: High-resolution spatio-temporal behavior of Parkfield seismicity[END_REF]earthquakes (Bakun et al., 2005;;[START_REF] Bakun | Earthquakes near Parkfield, California: Comparing the 1934 and 1966 Sequences[END_REF]. This behavior is also expected from results of dynamic modeling [START_REF] Das | An elliptical asperity in shear: fracture process and seismic radiation[END_REF] of rupture of an elliptical asperity, where the stress concentrations reach a maximum at the vertices of the semi-major axes of the ellipse.

The rupture terminates shortly after it impinges on the velocity-strengthening region.

This constitutes a stability barrier, as demonstrated by [START_REF] Kaneko | Towards inferring earthquake patterns from geodetic observations of interseismic coupling[END_REF]. However, recent studies [START_REF] Noda | Stable creeping fault segments can become destructive as a result of dynamic weakening[END_REF] show that some earthquake ruptures might also propagate through creeping areas that were considered to be barriers. The similarity of slip distribution of Parkfield earthquakes seems to suggest that the creeping zones stop rupture propagation at least for M6 and smaller earthquakes. The observation of logarithmically decaying afterslip sequences in the creeping section above the stuck patch following both the 1966 and 2004 earthquakes (e.g. [START_REF] Barbot | Postseismic deformation due to the Mw 6.0 2004 Parkfield earthquake: Stress-driven creep on a fault with spatially variable rateand-state friction parameters[END_REF][START_REF] Langbein | Coseismic and initial postseismic deformation from the 2004 Parkfield, California, earthquake, observed by global positioning system, electronic distance meter, creepmeters, and borehole strainmeters[END_REF][START_REF] Murray | Slip on the San Andreas fault at Parkfield, California, over two earthquake cycles, and the implications for seismic hazard[END_REF][START_REF] Scholz | Seismic and aseismic slip on the San Andreas Fault[END_REF] is consistent with the model of [START_REF] Marone | On the mechanics of earthquake afterslip[END_REF], in which a velocity strengthening region overlies a velocity weakening one.

Conclusions

We present a model based on geological surface observations and seismological observations at depth in which the three-dimensional structural complexities near MM and GH of the San Andreas Fault originate from the presence of a stuck patch surrounded by a creeping zone at Parkfield. The progressive helicoidal twist of the fault plane is due to a torque arising from the contrast in frictional strength between the stuck patch and the surrounding weak creeping area. The complex structures at each end of the helicoid result from the stress concentrations at these locations. They are characterized by complex reverse and normal faulting at the surface and diffuse seismicity at depth. These structural complexities are permanent features, playing an important role in stress concentrations and occurrence of M6 earthquakes at Parkfield. Faults are laterally segmented along their traces. In the early stages of their evolution, the fault segments are separated by geometrical discontinuities possibly formed by small secondary faults and fractures As the fault grows and lengthens through time, fault segments tend to become more and more connected and the lateral discontinuities in between smoothed out (e.g. Manighetti et al., 2015 and references therein). The main strand of the San Andreas fault at Parkfield is extraordinarily smooth, as a result of its large cumulative slip (Wesnousky, 1988); the 2-km fault offset at GH-CH being a notable exception. We find that this offset is a constructive edifice resulting in the strength contrast in the fault at this locality. This example of fault roughening with slip is in contrast to the generally assumed smoothing of fully locked faults at they become more mature.

ACKNOWLEDMENTS

This work is supported by the Brinson Foundation, and the National Science Foundation (NSF-EAR #1520680). We thank A. Lomax for his help with the NLLoc software. We California seismicity using cross-correlation and double-difference methods. J.

Fig

  Fig. S4 and S5 for detailed event distribution). It confirms that MM and GH-CH are

Figure 1 :

 1 FIGURES CAPTIONS

Figure 2 :

 2 Figure 2: a) Rotated map view of all surface fault traces from fig.1 (black lines) and the

Figure 3 :

 3 Figure 3: Along-strike cross-sections showing projected earthquake density (contours,

Figure 4 :

 4 Figure 4: (a) Geometry of the model domain and the fault plane. (b) Boundary

Figure 5 :

 5 Figure 5: (a) Perpendicular displacement of the fault plane toward the northeast and

Figure 6 :

 6 Figure 6: (a) Locations of the two vertical cross-sections corresponding to Middle

Figure 7 :

 7 Figure 7: Along-strike shear stress map on the fault plane at two stages of the

Figure 8 :

 8 Figure 8: Schematic map view (a) and cross section (b) showing the clockwise rotation

Figure 9 :

 9 Figure 9: Sketch showing a simplified view of the Parkfield fault section and along dip

thank the editor J.-P. Avouac, two anonymous reviewers and S. Barbot for their thorough comments that greatly improved the manuscript. Earthquake catalog data is available at http://ddrt.ldeo.columbia.edu. The DEM data (Lidar, SRTM3) have been provided by the OpenTopography Facility with support from the National Science Foundation under NSF Award Numbers 1226353 & 1225810. This paper is LDEO contribution 8319 and IPGP contribution 4043.


L Guironnet

P.-M Geffroy

N Tessier-Doyen

A Boulle

N Richet

T Chartier

P.-M Geffroy

Pierre-Marie Geffroy
email: pierre-marie.geffroy@unilim.fr

The surface roughness effect on electrochemical properties of La0.5Sr0.5Fe0.7Ga0.3O3-δ perovskite for oxygen transport membranes

Keywords: surface roughness, oxygen surface exchanges, oxygen semi-permeation, mixed electronic ionic conductors

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The oxygen transport membranes (OTM) are dense ceramic membranes usually manufactured using mixed ionic and electronic conductor materials. The electromotive force of the oxygen transport through the membrane is linked to the difference of oxygen partial pressures imposed between the two atmospheres on either side of the 2 membrane. However, the oxygen flux produced by the ceramic membranes is not enough high for an industrial development [START_REF] Sunsaro | Mixed ionicelectronic conducting (MIEC) ceramic-based membranes for oxygen separation[END_REF][START_REF] Geffroy | Rational selection of membrane materials of MIEC materials in energy production processes[END_REF]. Numerous studies reported in the literature aim to better understand the main mechanisms of oxygen transport through the membrane, in order to improve their performances of oxygen semi-permeation.

The majority of these studies reported that the oxygen flux though the membrane is governed by the kinetic of oxygen surface exchanges at the oxygen lean surface of the membrane. In this respect, a coating of ultra-divided particles or a porous coating on the membrane surface significantly improves oxygen flux through the membrane [START_REF] Hayamizu | Effects of surface modification on the oxygen permeation of Ba0.5Sr0.5Co0.8Fe0.2O3-δ[END_REF][START_REF] Guironnet | Improvement of oxygen flux through perovskite membranes using a coating of ultra-divided particles[END_REF][START_REF] Guironnet | La0.6Sr0.4Fe0.8Co0.2O3-δ electrophoretic coating for oxygen transport membranes[END_REF].

The effect of the roughness of the membrane surface on the oxygen semi-permeation performances is rarely studied. Tereaoka [START_REF] Kusaba | Surface effect on oxygen permeation through dense membrane of mixed-conductive LSCF perovskite-type oxide[END_REF] or Gurauskis et al. [START_REF] Gurauskis | La0.2Sr0.8Fe0.8Ta0.2O3δ based thin film membranes with surface modification for oxygen production[END_REF][START_REF] Ziesche | Permeation and oxygen exchange of Ln2Ni0.8Cu0.2O4-materials[END_REF][START_REF] Hoon | Oxygen Permeation and Stability of Ba0.5Sr0.5Co0.8Fe0.2O3-δ Membrane According to Trace Elements and Oxygen Partial Pressure in Synthetic Air[END_REF][START_REF] Bouwmeester | Importance of the surface kinetics as rate limiting step in oxygen permeation through mixed-conducting oxides[END_REF] show that the oxygen semi-permeation flux through La0.1Sr0.9Fe0.1Co0.9O3-δ and La0.2Sr0.8Fe0.8Ta0.2O3-δ perovskite membranes increases with the surface roughness of the membrane. Indeed, the surface quality of the membrane can lead to a variation of the oxygen semipermeation flux, if the oxygen flux through the membrane is governed by oxygen surface exchanges kinetics [START_REF] Kusaba | Surface effect on oxygen permeation through dense membrane of mixed-conductive LSCF perovskite-type oxide[END_REF][START_REF] Bouwmeester | Importance of the surface kinetics as rate limiting step in oxygen permeation through mixed-conducting oxides[END_REF]. Then, Wang [START_REF] Wei | Enhancement of oxygen permeation through U-shaped K2NiF4-type oxide hollow fiber membranes by surface modifications[END_REF] and Shaomin et al. [START_REF] Ning | Effect of enhanced oxygen reduction activity on oxygen permeation of La0.6Sr0.4Co0.2Fe0.8O3-δ membrane decorated by K2NiF4-type oxide[END_REF] reported that the surface treatment by acid etching of K2NiF4-type oxide membranes leads to improve the oxygen permeation performances and the oxygen surface exchange kinetics.

Besides, R.J. Carter et al. show than the apparent coefficient of oxygen surface exchange (k*), determined by the isotopic exchanges method, usually increases with the roughness of sample surface [START_REF] Chater | Development of a novel SIMS technique for oxygen self-diffusion and surface exchange coefficient measurements in oxides of high diffusivity[END_REF]. Then, the determination of k* obtained by the isotopic exchanges method suggests to control the surface finish of the sample. This work aims to evaluate the influence of the surface roughness on the oxygen flux through La0.5Sr0.5Fe0.7Ga0.3O3-δ perovskite-type membrane (LSFG5573) and, in particular on the kinetic coefficient of oxygen surface exchange of the membrane materials. The LSFG5573 material is of great interest as membrane material, with the advantages of good chemical stability [START_REF] Vivet | La(1-x)SrxFe(1-y)GayO3-δ perovskite membrane, oxygen semi-permeation, thermal expansion coefficient and chemical stability under reducing conditions[END_REF] in long duration tests, mechanical properties and high oxygen diffusion coefficient. Nevertheless, the main drawback of this material is a very low kinetic coefficient of oxygen surface exchange [START_REF] Ziesche | Permeation and oxygen exchange of Ln2Ni0.8Cu0.2O4-materials[END_REF].

Experimental part

Materials

The La0.5Sr0.5Fe0.7Ga0.3O3-δ (LSFG5573) perovskite-type powder is synthesized using a solid state reaction synthesis route. Pure oxides (La2O3 (99.99 %, Alfa Aesar), Fe2O3 (98 %, Alfa Aesar), Ga2O3 (99.999 %, Alfa Aesar)) are dehydrated at 100 °C. Then, dry oxides and carbonate (SrCO3 (99.9 %, Sigma-Aldrich)) precursors are weighed in the stoichiometric amount, mixed and attrition-milled at 600 rpm in ethanol media. The attrition milling is performed using 0.8 mm diameter zirconia balls. After 3 h, the mixture is separated from the zirconia balls with a 200 µm sieve, dried and calcined 5h at 1100°C. After calcination step, the perovskite phase is obtained, and the powder is once again attrition-milled in ethanol media during 3 hours to obtain a monomodal grain size distribution centered on 1 µm (laser granulometer -Horiba LA950). The chemical composition of the synthesized powder is determined by inductively coupled plasma atomic emission spectrometer (ICP-AES).

Elaboration of membrane samples

The membrane samples are elaborated by tape casting, as reported in previous works [START_REF] Geffroy | Identification of the ratedetermining step in oxygen transport through La(1-x)SrxFe(1-y)GayO3-δ perovskite membranes[END_REF][START_REF] Chartier | Tape casting[END_REF]. A LSFG5573 suspension is prepared from the synthesized powder, an organic solvent, a dispersant (phosphate ester), a binder (polymethacrylate) and a plasticizer (dibutyl phalate). The synthesized powder is first deagglomerated/dispersed for 1h in the organic solvent with the dispersant using ball-milling at 210 r.p.m. Then, the binder and the plasticizer are added and mixed to the slurry using ball-milling at 120 r.p.m for 12 h. The thickness of green tapes obtained by tape casting is 150 µm. The green tape is stamped into 30 mm diameter disks. Then, the disks are stacked and thermo-laminated under 50 MPa at 70 °C (in order to obtain 1 mm thick sintered disks). The green membranes are debinded and sintered at 1330 °C for 2h at a rate of 1 °C/min to obtain a relative density larger than 95% (Archimedes method).

Surface finishing of membrane samples

The surfaces of LSFG5573 membranes are polished with 3 various abrasive papers (SiC) in order to obtain different surface roughness. The "P120" membrane is roughly polished with a P120 grade-type abrasive paper, the "P600" membrane is polished with a softer abrasive paper and the "MP" membrane is mirror-polished. Each membranes were cleaned in an ethanol ultrasonic bath for ten minutes to eliminate the residual particles of silicon carbide. The roughness of the surfaces is measured by white light interferometry microscope (Fogale nanotech). The final thickness of membrane sample is adjusted during the polishing to be close to 1mm (+/-0.05 mm).

Structural characterizations of the membrane surface are carried out using grazing incidence X-ray diffraction (Bruker D8 Discover diffractometer, equipped with a copper source). A parabolic mirror coupled to a 2 reflections asymmetrical Germanium monochromator allows to select the Cu Kα1 wavelength and produces a parallel bean. The sample is positioned on a 5-axis goniometer. The detection is carried out using a LynxEyetype linear detector, covering an angular range of 2° with a resolution of 0.01°.

The oxygen semi-permeation performances of membrane samples

2.4.1

The setup for oxygen semi-permeation measurements Electrochemical performances, i.e. oxygen semi-permeation flux, and oxygen activities at the membrane surfaces are measured using a specific setup given Figure 1 and described in previous works [START_REF] Vivet | New route for high oxygen semi-permeation through surface-modified dense La1-xSrxFe1-yGayO3-δ pérovskite membranes[END_REF][START_REF] Geffroy | The effect of experimental factors on oxygen semi-permeation measurements[END_REF]. The membrane is sealed between two alumina tubes by gold rings [START_REF] Vivet | Influence of glass and gold sealants materials on oxygen permeation performances in La0.8Sr0.2Fe0.7Ga0.3O3-δ perovskite membranes[END_REF]. In the first chamber, recombined air is injected (100 mL/min, oxygen-rich atmosphere), in the second chamber, argon is injected (200 mL/min, oxygenlean atmosphere). The oxygen partial pressure gradient through the membrane or between the chamber 1 and 2, is the driving force for oxygen diffusion through the membrane. The chamber 3 prevents the eventual leakages between the chamber 1 and 2. The oxygen semi-permeation measurements are performed from 970 °C to 600°C

with an interval of 10 °C between each measure.

Figure 1. Oxygen semi-permeation setup and its microelectrodes system 2.4.2 Oxygen semi-permeation measurements

The oxygen flux through the membrane is evaluated using the following equation:

) (Eq. 1)
where J is the oxygen flux through the membrane in mol.m -2 .s -1 , D the argon flow in l.s -1 , V the molar volume of argon in l.mol -1 and S the effective membrane surface in m 2 . P and P ! are the oxygen partial pressure in the outlet argon flux and the oxygen partial pressure in the argon flux (in atm), respectively. They are evaluated from the electromotive force E(g) 2out, thanks to YSZ-oxygen sensors, following the Nernst law:

E # $% &' ln * +,- * . (Eq. 2)
where R is the universal gas constant, F the Faraday constant, T the temperature and PO2 is the oxygen partial pressure in the air (= 0.21 atm.).

The particularity of this apparatus is the possibility to determine the oxygen activity at the both surfaces of the membrane. The microelectrodes system, composed of a gold electrode and a zirconia tip, in contact with the surface measures two electromotive forces Es1 and Es2, respectively the variation of oxygen activities between the inlet gas and the oxygen-rich or lean surfaces of the membrane [START_REF] Fouletier | Electrochemical semipermeability and the Electrode Microsystem in Solid Oxide Electrolyte Cell[END_REF]. The gradients of di-oxygen activity on each surface of the membrane are evaluated from the Nernst law: (Eq. 4)

E /0 $% &' ln
where ∆μ 5

?!1@ μ 5

?!1@ μ A ?!1@ and ∆μ ?!1@ are the oxygen chemical potentials on oxygen-rich and oxygen-lean surfaces, respectively.

The difference between the oxygen chemical potential μ 5 ?!1@ and μ 5 BCD corresponds to the gradient of oxygen chemical activity in the bulk of the membrane, noted ∆μ E ?F .

Determination of oxygen surface exchange coefficients

ki and kdes are respectively the kinetic coefficients of oxygen surface exchange at oxygen-rich surface and at oxygen-lean surface of the membrane, expressed in cm.s -1 . They are described from the formalism recently reported by Bazan et al. [START_REF] Bazant | Theory of chemical kinetics and charge transfer based on nonequilibrium thermodynamics[END_REF] or Adler et al. [START_REF] Adler | Mechanisms and rate laws for oxygen exchange on mixed-conducting oxide surfaces[END_REF] and calculated from the oxygen semi-permeation flux and the gradient of oxygen chemical potential at the membrane surface, as following:

k B 2J IC K BCD Lexp P 0Q@ ∆R . 5 6789 $% S exp P Q@∆R . 5 6789 $% STU V (Eq. 5) k W!X 2J IC K ?!1@ Lexp P 0Q@ ∆R . 5 :;<= $% S exp P Q@∆R . 5 :;<= $% STU V (Eq. 6)
where n is a coefficient typically equal to 0.5 for MIEC materials and C K BCD and C K ?!1@ are the oxygen molar concentrations at the oxygen-rich and oxygen-lean surface, respectively. C is estimated by C

3 δ V ⁄ ,
where Vm is the molar volume estimated by XRD patterns and δ is the oxygen under-stoichiometry evaluated by thermogravimetric analysis under air.

The coefficient of oxygen diffusion, expressed in cm 2 s -1 , is calculated with the following equation:

D K &\] .
^_∆R .

`,:a (Eq. 7)

Identification of the limiting step of the oxygen transport through the membrane

Oxygen transport through the membrane is governed by bulk diffusion or/and kinetics of surface exchanges. In recent works, Valentin et al [START_REF] Valentin | Themo-chemico-mechanical modelling of mixed conductors in transient stages[END_REF][START_REF] Geffroy | Understanding and identifying the oxygen transport mechanisms through mixed-conductor membranes[END_REF] have introduced a new criterion called Bc, which is the ratio between the chemical potential gradient at the membrane surface and the one associated to the bulk diffusion. This criterion allows to identify the limiting step of oxygen transport through the membrane.

B C ∆R . 5

:;<= ∆R .

`,:a (Eq. 8)

If Bc is lower than 0.5, the oxygen flux is mainly governed by the oxygen diffusion through the membrane, and if Bc is higher than 1.5, the oxygen flux is mainly governed by oxygen surface exchanges between the gas and the membrane surfaces. Finally, in case of Bc is between 0.5 and 1.5, this corresponds to a mixed regime where the oxygen flux is partially governed by oxygen diffusion and surface exchanges. This suggests that the oxygen flux through the LSFG5573 membranes is mainly governed by the surface exchanges. The best oxygen semi-permeation flux is 3.9 10 -3 mol.m -².s -1 at 900°C, and it is obtained with "P120" membrane. This value of oxygen flux is close to the best performances reported in the literature with a dense membrane thickness of 1 mm. For instance, the best oxygen semi-permeation fluxes reported in the literature with Ba1-xSrxFe1-yCoyO3-δ perovskite membranes are from 1 to 20 10 -3 mol.m -2 .s -1 at 900 °C [START_REF] Ge | Oxygen selective membranes based on B-site cation-deficient (Ba0.5Sr0.5)(Co0.8Fe0.2)yO3-δ perovskite with improved operational stability[END_REF]26].

Results and discussion

The dependence of the surface quality on the oxygen semi-permeation

Figure 3 shows the variation of Bc lean criterion in relation with the temperature and the surface roughness of the membrane. Bc lean is higher to 1.5 for "As sintered" and "MP" membranes, corresponding to an oxygen flux governed by surface exchanges at the oxygen lean side of the membrane, as expected. In case of P600 and P120 membranes, the variation of Bc lean with the temperature shows that the oxygen flux is partially governed by oxygen diffusion at high temperature (> 900 °C), with a Bc lean close to 1. In opposite case, the oxygen flux is mainly governed by surface exchange at lower temperature (< 800 °C), with a Bc lean > 1.5. In consequence, the Arrhenius plots of the oxygen semi-permeation flux through the "P600" and "P120" membranes show a variation of activation energy (figure 2b), which is likely linked to a change of the limiting step of oxygen transport through the membrane with the temperature, as reported by the evolution of Bc lean in figure 3. Besides, the activation energy is decreasing with the surface roughness of the LSFG5573 membrane, from 130 kJ.mol -1 for "As sintered" membrane to 190kJ.mol -1 for "PM" membrane at high temperature. This suggests that there is an evolution of the surface exchange mechanism as function of the surface state of the membrane. This leads to understand the influence of the surface physico-chemistry on the kinetics of oxygen surface exchanges, as reported in the next section.

Understanding of the surface physico-chemistry on the kinetics of oxygen surface exchanges

The roughness of the membrane surfaces (Ra) is evaluated by interferometric microscopy. It is possible to design 3D and 2D cartographies of the surface topographies (Fig. 4). The software (Mountains Software) allows to calculate the real surface developed by the roughness from topographic data of the membrane surface. The surface developed by the roughness is determined here from the surface ratio (Sd/Sapp) between the surface developed by the roughness (Sd) and the apparent surface of the membrane (Sapp) (Table 1). The variation of this ratios (Sd/Sapp) remains low for the four surfaces characterized, i.e. between 1.1 and 1.4. This means that the real surface developed by the roughness is still close to the apparent surface of the membrane in this study.

Tereoka et al. reported significantly higher ratios in relation with the roughness of the membrane surface, i.e. between 1.44 and 3.19 [START_REF] Kusaba | Surface effect on oxygen permeation through dense membrane of mixed-conductive LSCF perovskite-type oxide[END_REF].

Table 1 gives the average surface roughness in relation with the abrasive paper grade used for the polishing of the membrane surface. Logically, a surface polished with a P120 paper has higher roughness (Ra= 0.86 µm) than one polished with a P600 paper (Ra= 0.54 µm). The average roughness of the mirror-polished surface membrane is very low (Ra= 0.15 µm) in comparison with other surfaces. The roughness of the membrane surface increases as following: MP < P600 < As sintered < P120. Figure 5 shows the XRD patterns (obtained at grazing angle of 5°) recorded from the membranes with the different polishing grades. Secondary phases, such as SrLaGa3O7, are very difficult to totally eliminate by synthesis of LSFG5573 by solid state reaction [START_REF] Geffroy | Identification of the ratedetermining step in oxygen transport through La(1-x)SrxFe(1-y)GayO3-δ perovskite membranes[END_REF][START_REF] Kharton | Perovskite-like system (Sr,La)(Fe,Ga)O3-δ: Structure and ionic transport under oxidizing conditions[END_REF] and this phase is observed in the as-sintered and for the first polishing state. An important outcome of these experiments is that they evidence a segregation of the SrLaGa3O7 secondary phase at the surface of LSFG5573 membrane during the sintering. Indeed, after intense polishing of the membrane surface, the SrLaGa3O7 secondary phase is no more detected on XRD patterns, because the concentration of the SrLaGa3O7 secondary phase in bulk is too low to be detected on XRD patterns.

This secondary phase is only detected at the vicinity of the surface by XRD analysis at grazing angle after sintering. Then, the surface polishing of the membrane leads to remove (totally or partially) this secondary phase at the membrane surface. From the relative intensities of the 110 peak of LSFG5573 (2θ = 32.46°) and 211 peak of SrLaGa3O7 (2θ = 29.91°), and using the integrated intensity equation [START_REF] Warren | X-ray diffraction[END_REF], it is possible to derive the volumic ratio of those phase. The results are shown in Table 2. Figure 6 shows the effect of the roughness/secondary phase on the oxygen semi-permeation fluxes through the LSFG5573 membranes at 900 °C. The oxygen semi-permeation fluxes obtained for "as sintered" membrane and "P600" membrane are of the same order of magnitude, i.e. 8.5 and 12 10 -4 mol.m -2 .s -1 respectively at 900 °C. It is slightly lower, i.e. 2.2 10 -4 mol.m -2 .s -1 , for the membrane with a mirror-polished surface. The larger flux is obtained for the "P120" membrane with 3.9 10 -3 mol.m -2 .s -1 at 900 °C. The variation of the oxygen semipermeation flux with the modification of the surface characteristics suggests that the oxygen surface exchange at the membrane surface is the limiting step of oxygen transport through the membrane, as discussed in the previous section.

Two contributions of the surface on the kinetics of oxygen surface exchanges can been considered to explain the evolution of oxygen flux: i) the roughness and ii) the quantity of secondary has a large effect. Considering the first possible contribution, the variation of the membrane surface with the roughness is very low (< 1.4) in comparison with the variation of oxygen flux (an order of magnitude higher). This suggests that the variation of the exchange surface due to the roughness, in the range tested, has a low effect on the kinetics of oxygen surface exchange.

Then, the variation of oxygen flux with the different surfaces may likely be attributed mainly to the segregation of SrLaGa3O7 secondary phase during sintering. SrLaGa3O7 is reported as a dielectric material in the literature [START_REF] Warren | X-ray diffraction[END_REF][START_REF] Xu | Oxygen interstitials and vacancies in LaSrGa3O7-based melilites[END_REF][START_REF] Huang | Superior perovskite oxide-ion conductor; strontium-and magnesium-doped LaGaO3[END_REF], which leads to decrease the mobility of electrons or oxygen anions at the surface of mixed electronic and ionic conductors. membrane is far from the equilibrium state in oxygen semi-permeation conditions. In this respect, the surface reaction can be associated to incorporation or desorption on oxygen rich and lean sides, as reported in reaction 1 and 2: The coefficient of oxygen incorporation at the oxygen-rich side (ki) of the membrane exponentially increases with the roughness of the membrane surface. There is not apparent effect of the SrLaGa3O7 secondary phase on the kinetics of incorporation of oxygen.

O 2 g + V O ⦁⦁ F g → O O X + 2 ℎ ⦁
However, the oxygen flux through the LSFG5573 membranes is predominantly governed by oxygen surface exchanges at the oxygen-lean side of the membrane attested by a ki coefficient significantly larger than kdes coefficient for the same roughness values (Ra), 1.6 10 -6 and 3.6 10 -7 cm.s -1 respectively, for P600 membrane at 900 °C (Fig. 7).

The presence of SrLaGa3O7 secondary phase on the membrane surface (only for the sample noted "as sintered") seems to have a large effect on the coefficients of oxygen desorption at the oxygen lean side of the membrane.

The kdes coefficient of the "as sintered" membrane is lower than the kdes coefficient of the "P600" membrane, 4.1 10 -7 and 1.7 10 -6 cm.s -1 respectively, although, the roughness (Ra) of the "as sintered" membrane surface is higher than the one of the "P600" membrane surface. Besides, the SrLaGa3O7 secondary phase is a dielectric material [START_REF] Warren | X-ray diffraction[END_REF][START_REF] Xu | Oxygen interstitials and vacancies in LaSrGa3O7-based melilites[END_REF][START_REF] Huang | Superior perovskite oxide-ion conductor; strontium-and magnesium-doped LaGaO3[END_REF], which leads to decrease the mobility of electrons or oxygen species at the surface of mixed conductors, and then the kinetics of oxygen desorption at the membrane surface.

Finally, the large effect of the surface roughness on the kinetic of oxygen desorption or incorporation is not well understood here. For instance, XPS analysis at the membrane surface shows not significant impurities or chemical segregations for the different surfaces reported in this work. The main assumption is that a large roughness could be favorable to increase the concentration of oxygen vacancies at the membrane surface. The oxygen vacancy concentration could be locally modified by the local surface radius at the interface solid/gas, as reported by Laplace equation in the literature on the sintering mechanism [32]. This local variation of oxygen vacancy concentration at the material surface could strongly increases the kinetics of oxygen surface or the apparent coefficients of oxygen exchanges at the membrane surfaces.

Conclusions

The surface roughness of LSFG5573 perovskite membrane has a large effect on their oxygen semi-permeation performances. The oxygen flux through LSFG5573 perovskite membrane is mainly governed by the kinetics of oxygen desorption at the oxygen lean side of the membrane.

Besides, the presence of SrLaGa3O7 secondary phase at the LSFG5573 perovskite membrane largely decreases the coefficient of oxygen desorption at the oxygen lean sides. This SrLaGa3O7 secondary phase segregates at the membrane surface during sintering, and the surface polishing of the membrane leads to remove (totally or partially) this secondary phase at the membrane surface.

Then, a large roughness at the membrane surface without secondary phase leads to an increase of the kinetics of oxygen surface exchanges, and indirectly the oxygen flux through the membrane. By this way, the best oxygen fluxes obtained in this work, is 3.9 10 -3 mol.m -2 .s -1 , and is not far from the best performances reported in the literature with a dense 1 mm in thick membrane.

These results require to define correctly the roughness and the physico-chemistry of the membrane surface, to evaluate the oxygen semi-permeation of the membrane materials. Finally, this work shows that it is possible to largely improve the electrochemical properties of membranes by a control of the surface physico chemistry. This could lead to a promising route for the development of new design for cathode materials for solid oxide fuel cell or for membranes for oxygen separation applications.

Table captions

Table 1. Roughness at the membrane surfaces obtained by interferometric microscopy Table 2. Qualitative estimation of SrLaGa3O7 type secondary phase as function of the surface roughness

  and μ A ?!1@ are the oxygen chemical potentials in air and argon. μ 5BCD and μ 5

Figure 2

 2 Figure 2 shows, a) the oxygen semi-permeation flux as function of the temperature for LSFG5573 membranes and, b) the Arrhenius plots of oxygen semi-permeation for LSFG5573 membranes with different surface roughness.

Figure 2 .

 2 Figure 2. The oxygen semi-permeation fluxes through LSFG5573 membranes a) function of the temperature and b) function of the surface roughness at 900 °C

Figure 3 .

 3 Figure 3. Variation of Bc lean with temperature for different surface roughness of membranes

Figure 4 .

 4 Figure 4. Cartographies of the membrane surfaces obtained by white light interferometry, A) As sintered surface, B) P120, C) P600 and D) MP.

Figure 5 .

 5 Figure 5. XRD patterns at grazing angle of LSFG5573 membranes as function of the surface roughness

Figure 6 .Figure 7

 67 Figure 6. Effect of surface roughness on the oxygen flux through the membrane at 900 °C Figure 7 shows the coefficients of oxygen surface exchanges at 900 °C at the oxygen rich and lean sides, in relation with the roughness (Ra). ki and kdes coefficients can be attributed to the surface exchange coefficients of oxygen incorporation and oxygen desorption on each side of the membrane, because the surface of LSFG5573

reaction 1 :Figure 7 .

 17 Figure 7. Effect of the roughness on the coefficients, a) of oxygen incorporation at the oxygen-rich side (ki) and b) of oxygen desorption at the oxygen lean side (kdes) at 900 °C.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure 1. Oxygen semi-permeation setup and its microelectrodes systemFigure 2. The oxygen semi-permeation fluxes through LSFG5573 membranes a) function of the temperature and b) function of the surface roughness at 900 °C Figure 3. Variation of Bc lean with temperature for different surface roughness of membranes Figure 4. Cartographies of the membrane surfaces obtained by white light interferometry, A) As sintered surface, B) P120, C) P600 and D) MP.

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. XRD patterns at grazing angle of LSFG5573 membranes as function of the surface roughness Figure 6. Effect of surface roughness on the oxygen flux through the membrane at 900 °C Figure 7. Effect of the roughness on the coefficients, a) of oxygen incorporation at the oxygen-rich side (ki) and b) of oxygen desorption at the oxygen lean side (kdes) at 900 °C.

Table 1 .

 1 Roughness at the membrane surfaces obtained by interferometric microscopy

	Polishing	P120	As sintered	P600	MP
	Ra (µm)	0.86	0.54	0.38	0.15
	Sd/Sapp	1.4	1.2	1.1	1.1

Acknowledgements

The authors wish to express their gratitude to Air Liquide Company for financial and technical support and the ADEME (French Environment and Energy Management Agency) for financial support of this study.