Louis De Barros

Bastien Dupuy

Gareth S O'brien

Jean Virieux

Stéphane Garambois

Using a Poroelastic Theory to Reconstruct Subsurface Properties: Numerical Investigation

Introduction

The quantitative imaging of the Earth subsurface is a major challenge in geophysics. In oil and gas exploration and production, aquifer management and other applications such as the underground storage of CO 2 , seismic imaging techniques are implemented to provide as much information as possible on fluid-filled reservoir rocks. Biot theory [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range, II. Higher frequency range[END_REF] and its extensions provide a convenient framework to connect the various parameters characterizing a porous medium to the wave properties, namely, their amplitudes, velocities and frequency contents. The poroelastic model involves more parameters than the elastodynamic theory, but on the other hand, the wave attenuation and dispersion characteristics at the macroscopic scale are determined by the intrinsic properties of the medium without having to resort to empirical relationships. Attenuation mechanisms at microscopic and mesoscopic scales, which are not considered in the original Biot theory, can be introduced into alternative poroelastic theories (see e.g. [START_REF] Pride | Seismic attenuation due to wave-induced flow[END_REF].

The inverse problem, that is, the retrieval of poroelastic parameters from the seismic waveforms, is much more challenging. Porosity, permeability and fluid saturation are the most important parameters for reservoir engineers. The estimation of poroelastic properties of reservoir rocks from seismic waves is however still in its infancy. The classical way of estimating these is to first solve the elastic problem and then interpret the velocities in terms of poroelastic parameters by using deterministic or stochastic rock physics modelling. However, unlike Full Waveform Inversion (FWI), these methods do not make full use of the seismograms.

In the poroelastic case, eight model parameters are required to describe the medium, compared with only one or two in the acoustic case, and three in the elastic case if wave attenuation is not considered. The advantages of using a poroelastic theory in FWI are (1) to directly relate seismic wave characteristics to porous media properties; (2) to use information that cannot be described by viscoelasticity or elasticity with the [START_REF] Gassmann | Über die elastizität poröser medien[END_REF] formulae and ( 3) to open the possibility to use fluid displacement and force to determine permeability and fluid properties.

As an example of geological target, figure 1 presents the data recorded by a seismic survey on a seashore in the South of France. As water is pumped inland, saltwater is intruding into the coastal aquifer. This can affect the ground water and lead to severe problems with water supplies in the area. The monitoring of this phenomenon requires knowledge of the soil characteristics, including the permeability and porosity, and the properties of the fluid. A simple elastic approach cannot fully solve this problem, however the poroelastic theory may offer an alternative solution. In this example, the medium comprises alternating layers of sand, silt and clay with varying levels of compaction and a wide range of porosity and permeability. This layering produces strong reflected waves as shown in figure 1. The aim of the paper is to investigate how a poroelastic theory can be used to monitor water flow, and identify preferential pathways, using reflected waves. However, at this early stage, this chapter will only focus on numerical tests. Example of data for which a poroelastic-based interpretation may be useful. Data are recorded on a seashore in the South of France. Source is a hammer shot, and the 24 receivers are equally spaced between 10 and 55 m from the source. The medium is very soft and water saturated, with a direct P wave velocity of c. 1600m/s and a S wave velocity lower than 200 m/s. We will investigate the gain and the feasibility of using a poroelastic approach, rather than the classical elastic one, in full waveform methods. The forward modelling is solved using different algorithms: a reflectivity approach, a 3D finite difference scheme and a 2D discontinuous Galerkin method. The comparison of synthetic data computed in the elastic and poroelastic cases shows that poroelastic modelling leads to some typical patterns that cannot be explained by elastic theory. This proves that the use of poroelastic theories may bring more insight to the model reconstruction, particularly, in relation to the fluid properties. Moreover, mesoscopic attenuation can be introduced in the poroelastic laws for double porosity medium, adding extra changes in the waveforms. This demonstrates the utility of using such theories to correctly reproduce measured seismic data.

Analytical formulae are then derived to compute the first-order effects produced by plane inhomogeneities on the point source seismic response of a fluid-filled stratified porous medium. The derivation is achieved by a perturbation analysis of the poroelastic wave equations in the plane-wave domain using the Born approximation. The sensitivity of the wavefields to the different model parameters can be investigated: the porosity, consolidation parameter, solid density, and mineral shear modulus emerge as the most sensitive parameters in the forward and inverse modelling problems. However, the amplitude-versus-angle response of a thin layer shows strong coupling effects between several model parameters.

The inverse problem is then tackled using a generalized least-squares, quasi-Newton approach to determine the parameters of the porous medium. Simple models consisting of plane-layered, fluid-saturated and poro-elastic media are considered to demonstrate the concept and evaluate the performance of such a full waveform inversion scheme. Numerical experiments show that, when applied to synthetic data, the inversion procedure can accurately reconstruct the vertical distribution of a single model parameter, if all other parameters are perfectly known. However, the coupling between some of the model parameters does not permit the reconstruction of several model parameters at the same time. To get around this problem, we consider composite parameters defined from the original model properties and from a priori information, such as the fluid saturation rate or the lithology, to reduce the number of unknowns. We then apply this inversion algorithm to time-lapse surveys carried out for fluid substitution problems, such as CO 2 injection, since in this case only a few parameters may vary as a function of time. A two-step differential inversion approach allows the reconstruction of the fluid saturation in reservoir layers, even though the medium properties are mainly unknown.

Wave propagation in stratified porous media

The governing equations for the poroelastodynamic theory were first derived by [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range, II. Higher frequency range[END_REF], and are thus often referred to as "Biot's theory". The main hypothesis behind these equations is that the seismic wavelengths are longer than the pore size; the medium can then be described by homogeneised laws. Poroelastic theories have since been derived and improved by many authors (e.g. [START_REF] Auriault | Dynamics of porous saturated media, checking of the generalized law of Darcy[END_REF][START_REF] Geertsma | Some aspects of elastic wave propagation in fluid-saturated porous solid[END_REF][START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF].

Governing equations

Assuming a e -iωt dependence, [START_REF] Pride | Deriving the equations of motion for porous isotropic media[END_REF] rewrote [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range, II. Higher frequency range[END_REF] equations of poro-elasticity in the frequency domain in the form

[ (K U + G/3) ∇∇ + (G∇ 2 + ω 2 ρ) I ] . u + [ C∇∇ + ω 2 ρ f I ] . w = 0 ( 1 ) [ C∇∇ + ω 2 ρ f I ] . u + [ M∇∇ + ω 2 ρI ] . w = 0 ,
where u and w respectively denote the average solid displacement and the relative fluid-to-solid displacement, ω is the angular frequency, I the identity tensor, ∇∇ the gradient of the divergence operator and ∇ 2 the Laplacian operator. The other quantities appearing in equations (1) are properties of the medium. The bulk density of the porous medium ρ is related to the fluid density ρ f , solid density ρ s and porosity φ:

ρ = (1 -φ)ρ s + φρ f . ( 2
)
K U is the undrained bulk modulus and G is the shear modulus. M (fluid storage coefficient) and C (C-modulus) are mechanical parameters. In the quasi-static limit, at low frequencies, these parameters are real, frequency-independent and can be expressed in terms of the drained bulk modulus K D , porosity φ, mineral bulk modulus K s and fluid bulk modulus K f [START_REF] Gassmann | Über die elastizität poröser medien[END_REF]:

K U = φK D + 1 -(1 + φ) K D K s K f φ(1 + Δ) , C = 1 - K D K s K f φ(1 + Δ) , M = K f φ(1 + Δ) (3) with Δ = 1 -φ φ K f K s 1 - K D (1 -φ)K s .
It is also possible to link the frame properties K D and G to the porosity and constitutive mineral properties [START_REF] Korringa | Self-consistent imbedding and the ellipsoidal model for porous rocks[END_REF][START_REF] Pride | Relationships between seismic and hydrological properties, in Hydrogeophysics[END_REF]:

K D = K s 1 -φ 1 + c s φ and G = G s 1 -φ 1 + 3c s φ/2 , ( 4
)
where G s is the shear modulus of the grains. The consolidation parameter c s appearing in these expressions is not necessarily the same for K D and G [START_REF] Korringa | Self-consistent imbedding and the ellipsoidal model for porous rocks[END_REF]. However, to minimize the number of model parameters, and following the recommendation of [START_REF] Pride | Relationships between seismic and hydrological properties, in Hydrogeophysics[END_REF], we consider only one consolidation parameter to describe the frame properties. c s typically varies between 2 to 20 in a consolidated medium, but can be much greater than 20 in a soft soil.

Finally, the wave attenuation is explained by a generalized Darcy's law which uses a complex, frequency-dependent dynamic permeability k(ω) defined via the relationship [START_REF] Johnson | Probing porous media with first and second sound. I. dynamic permeability[END_REF]:

ρ = i η ω k(ω) with k(ω) = k 0 / 1 -i 4 n J ω ω c -i ω ω c . ( 5 )
In equation ( 5), η is the viscosity of the fluid and k 0 the hydraulic permeability. Parameter n J is considered constant and equal to 8 to simplify the equations.

The relaxation frequency ω c = η/(ρ f Fk 0 ), with F the electrical formation factor, separates the low frequency regime where viscous losses are dominant from the high frequency regime where inertial effects prevail. We refer the reader to the work of [START_REF] Pride | Relationships between seismic and hydrological properties, in Hydrogeophysics[END_REF] for more information on the parameters used in this study.

The solution of equation (1) leads to classical fast P-and S-waves, and to an additional slow P-wave (often called Biot wave) . The fast P-wave has fluid and solid motion in phase, while the Biot wave has out-of-phase motions. At low frequency, the Biot wave has a diffusive pattern and can be seen as a fluid pressure diffusion wave. At high frequency, the inertial effects are predominant. This wave becomes propagative and can be seen in data, giving an experimental justification to the dynamic poroelasticity theory [START_REF] Plona | Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies[END_REF].

Mesoscopic attenuation and more complex theories

Although the slow P-wave does not appear on the seismograms at low frequency, it plays an important role in the attenuation process, as it produces loss of energy by wave-induced fluid-flow. However, the attenuation as described in the Biot theory is not strong enough to model the attenuation in geological media, especially at low (i.e. seismic) frequencies.

Attenuation processes can actually be separated into 3 different spatial scales, namely microscopic, mesoscopic and macroscopic [START_REF] Pride | Seismic attenuation due to wave-induced flow[END_REF]. Within this classification, the Biot mechanism of attenuation takes place at macroscopic scale (on the order of the seismic wavelength). The microscopic attenuation is due to mechanisms that occur at the grain size, such as the squirt flow mechanism [START_REF] Mavko | Estimating grain-scale fluid effects on velocity dispersion in rocks[END_REF]. This mechanism leads to wave attenuation mainly at high frequencies. The attenuation mechanism that prevails at low frequency comes from the mesoscopic scale [START_REF] Pride | Seismic attenuation due to wave-induced flow[END_REF], and it is due to fluid flow that occurs at boundaries between any medium heterogeneities whose sizes are between the grain sizes and the seismic wavelengths. This is particularly true for layered media [START_REF] Gurevich | Seismic attenuation in finely layered porous rocks: Effects of fluid flow and scattering[END_REF][START_REF] Pride | Biot slow-wave effects in stratified rock[END_REF] or when the medium contains 1) inclusions of different materials such as composite medium or double porosity medium [START_REF] Berryman | Elastic wave propagation and attenuation in a double-porosity dual permeability medium[END_REF][START_REF] Pride | Seismic attenuation due to wave-induced flow[END_REF][START_REF] Santos | On the static and dynamic behavior of fluid saturated composite porous solids: A homogenization approach[END_REF], or 2) different fluids [START_REF] Santos | A model for wave propagation in a porous medium saturated by a two-phase fluid[END_REF] or patches of different saturation [START_REF] Johnson | Theory of frequency dependent acoustics in patchy-saturated porous media[END_REF].

Double porosity medium (DP in the text) refers to a porous medium which contains inclusions with different porosity and permeability [START_REF] Pride | Seismic attenuation due to wave-induced flow[END_REF][START_REF] Santos | On the static and dynamic behavior of fluid saturated composite porous solids: A homogenization approach[END_REF]. Assuming that the most compressible phase (patches, phase 2) is embedded into the least compressible one (host rock, phase1), the fluid flow inside the phase 2 can be eliminated from the homogenised equations. This assumption allowed Pride & Berryman (2003a[START_REF] Pride | Linear dynamics of double-porosity dual-permeability materials, ii. fluid transport equations[END_REF] to write the DP equations under the form of the classical Biot theory (eq. 1). This involves the use of complex frequency dependent moduli K U , C and M. Particularly, these parameters are functions of the respective volume of each phase and of the size of the patches (a denotes here the average radius of the inclusions in the host rock). Finally, as the patches are assumed to be spherical, the shear modulus of the medium is still real and not frequency dependent, and can be approximated by the geometrical mean of the modulus inside each phase. For the derivation and the detailed expressions of the parameters, please refer to the work of [START_REF] Pride | Seismic attenuation due to wave-induced flow[END_REF] and [START_REF] Pride | Relationships between seismic and hydrological properties, in Hydrogeophysics[END_REF].

To show how the mesoscopic attenuation due to DP media impacts the seismic properties, we look for the changes in the P-wave velocity and attenuation. The medium is composed of little patches of high permeable and high porosity in a less permeable host rock. Following the work of [START_REF] Liu | Transient solution for poro-viscoacoustic wave propagation in double porosity media and its limitations[END_REF], we consider a sandstone with 3% sand inclusions. The complex moduli K U , C and M are computed using the DP effective theory of [START_REF] Pride | Seismic attenuation due to wave-induced flow[END_REF], leading to the P-wave velocity and attenuation with respect to the frequency. The results are compared to the seismic properties of each single phase and to the response using an average single porosity medium, where moduli are computed by geometrical averages. Figure 2 shows the P-wave velocity and attenuation (via the inverse of the quality factor) for the double porosity medium (for the inclusion radius a equal to 1, 5 and 10 cm), for the average single porosity medium and for both single phase media.

The P-wave velocity is much more dispersive for the double porosity medium than for the equivalent single porosity medium. At high frequency, it is much higher in the double porosity medium than in the equivalent single porosity medium. It shows two main changes: and10 8 Hz). Blue lines are for the double porosity theory for inclusion sizes a = 1 cm (continuous lines), a = 5 cm (dotted-dashed lines) and a = 10 cm (dashed lines). Poroelastic responses of the individual phases are the black lines. The results for the average single porosity theory, computed using weighted average of each phase properties, are given by the red lines. 1) at the relaxation frequency of the host rock medium, which are at the transition between the diffusive and the propagative regime of the Biot wave, and 2) at low frequencies, around the relaxation frequency of the patches, with a strong dependance on a. It is worth noting that the size of inclusions, a, has an strong effect on the low frequency behaviour. We observe similar patterns for the P-wave attenuation (inverse of the quality factor). In the double porosity medium, attenuation shows two main peaks, associated with the two phases, while the equivalent single poroelastic medium produces only one single peak at high frequency. At low frequency (seismic frequencies), the attenuation is very high for the double porosity medium, leading to quality factors that are in good agreement with quality factors in geological materials. This pattern strongly depends on the size of the inclusions. On the other hand, attenuation produced by single phase medium is too low to be realistic. This means that the fluid flow at the boundaries between heterogeneities plays a fundamental role in the attenuation process, that cannot be neglected. These P-wave characteristics will have a strong influence on the seismic waveforms (see section 3.2).

Using a poroelastic theory is much more complex than an equivalent visco-elastic theory. However, modelling seismic waves with poroelastic theories take into account the attenuation induced by fluid equilibration at layer interfaces or heterogeneity boundaries, whereas a viscoelastic approach neglects this attenuation process. As shown by [START_REF] Pride | Seismic attenuation due to wave-induced flow[END_REF], this is the most important attenuation process at low frequency. As the shallow subsurface has strong lateral and vertical heterogeneities, one should solve the full poroelastic theory to deal with attenuation.

Numerical modelling of seismic waves in porous media

Forward modelling solution

The model properties m, which are the material parameters introduced in the previous section, are nonlinearly related to the seismic data d via an operator f , i.e., d = f (m). The forward problem has been solved by many authors, using different methods. Analytical solutions have been derived for a homogeneous medium [START_REF] Boutin | Green functions and associated sources in infinite and stratified poroelastic media[END_REF][START_REF] Philippacopoulos | Buried point source in a poroelastic half-space[END_REF]. The response of porous layered medium has been computed using reflectivity methods in the frequency-wavenumber domain, such as the [START_REF] Kennett | Seismic Wave Propagation in Stratified Media[END_REF] approach [START_REF] De Barros | Perturbations of the seismic reflectivity of a fluid-saturated depth-dependent poroelastic medium[END_REF][START_REF] Pride | Biot slow-wave effects in stratified rock[END_REF]. This method was also used to solve the coupling between seismic and electromagnetic waves [START_REF] Garambois | Full waveform numerical simulations of seismoelectromagnetic wave conversions in fluid-saturated stratified porous media[END_REF][START_REF] Haartsen | Electroseismic waves from point sources in layered media[END_REF]. The poroelastic equations have been solved in 2D and 3D cases, mainly using finite difference schemes [START_REF] Carcione | Viscoelastic effective rheologies for modelling wave propagation in porous media[END_REF][START_REF] Dai | Wave propagation in heteregeneous porous media: A velocity-stress, finite-difference method[END_REF][START_REF] Masson | Finite-difference modeling of Biot's poroelastic equations across all frequencies[END_REF][START_REF] O'brien | 3D rotated and standard staggered finite-difference solutions to Biot's poroelastic wave equations: Stability condition and dispersion analysis[END_REF] in the time-space domain. For discretisation issues, the equations (1) should be decomposed in propagative and diffusive parts, which have to be solved independently [START_REF] Carcione | Viscoelastic effective rheologies for modelling wave propagation in porous media[END_REF]. Other time domain numerical schemes have been used, such as finite elements [START_REF] Morency | Spectral-element simulations of wave propagation in porous media[END_REF] or finite volume (de la [START_REF] De La Puente | Discontinuous galerkin methods for wave propagation in poroelastic media[END_REF]. Finally, [START_REF] Dupuy | Wave propagation in heterogeneous porous media formulated in the frequency-space domain using a discontinuous galerkin method[END_REF] solved this problem in the frequency domain using a discontinuous Galerkin approach. For a complete and precise review of the numerical modelling used to solve the poroelastic problem, we refer the reader to [START_REF] Carcione | Computational poroelasticity -a review[END_REF].

In this paper, we will use three different techniques to illustrate our points:

• a 3D Finite Difference scheme (FD; O'Brien, 2010): The solutions of Biot's equations are obtained by fourth-order in space and second-order in time staggered-grid and rotated-staggered-grid methods. Stability of the methods and accuracy of the solutions have been carefully checked in the low-frequency domain. • a reflectivity approach (SKB; De [START_REF] De Barros | Perturbations of the seismic reflectivity of a fluid-saturated depth-dependent poroelastic medium[END_REF][START_REF] De Barros | Full waveform inversion of seismic waves reflected in a stratified porous medium[END_REF]: The 3D solution is obtained in the frequency-wavenumber domain for horizontally layered media by using the generalized reflection and transmission method of [START_REF] Kennett | Seismic Wave Propagation in Stratified Media[END_REF]. The synthetic seismograms are finally transformed into the time-distance domain by using the 3D axisymmetric discrete wavenumber integration technique of [START_REF] Bouchon | A simple method to calculate Green's functions for elastic layered media[END_REF]. • a Discontinuous Galerkin Method (DGM; [START_REF] Dupuy | Wave propagation in heterogeneous porous media formulated in the frequency-space domain using a discontinuous galerkin method[END_REF]: For 2 dimension medium, the discrete linear system for the Biot theory has been deduced in the frequency domain for a discontinuous finite-element method, known as the nodal discontinuous Galerkin method. Solving this system in the frequency domain allows accurate modelling of the wave propagation for all frequencies.

The last two approaches are in the frequency domain, which has several advantages: 1) there is no need to decompose the problem into diffusive and propagative parts; 2) all frequencies, i.e., in the low-and high-frequency regimes, can be accuarately treated; 3) solving more complex theories, such as double porosity or poroviscoelastic theories is straighforward and does not require any modifications of the solver; and 4) frequency domain has been shown to be the most efficient way to solve the Full waveform inverse problem, as the solution has to be calculated only for a few frequencies [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF].

Seismic waveforms in poroelastic medium

As already stated, the main improvement of using poroelastic versus elastic theories is in the description of the attenuation from intrinsic medium parameters. Figure 3 gives an example of poroelastic and elastic data computed by equivalent finite difference codes (FD, O'Brien, 2010). The medium is a complex 200m thick reservoir embedded in a homogeneous half-space. The reservoir, modified from [START_REF] Manzocchi | Sensitivity of the impact of geological uncertainty on production from faulted and unfaulted shallow marine oil reservoirs -objectives and methods[END_REF] is composed of 7 different facies, with different mineral and frame properties (see Fig. 3,top). In order to mimic a time lapse survey for CO 2 geological storage in a saline aquifer, a baseline is first computed for fully brine saturated medium. CO 2 is then injected in the center of the reservoir and spread out according to the permeability, leading to areas containing gas of roughly 500 m and 2000 m diameter.

Figure 3 presents the differential data (data with CO 2 in the reservoir minus baseline) for both gas extensions. The same models are run using equivalent elastic properties (computed through the Gassmann formulation). As the velocities are equal, arrival times of the elastic and poroelastic waves are the same. However, changes appear in amplitude, mainly in the multiple reflected waves and coda. The amplitude differences are up to 40% of the data. It stresses the importance of using attenuation in the forward modelling and inversion processes. Even if the Biot theory cannot explain the full range of attenuation, it leads to significant changes in the seismic waveforms.

Poroelastic attenuation is due to the fluid movement. In the poroelastic theories, evidence for this lies in the relative fluid-to-solid motion and in the existence of the slow P-wave. Figure 4 reproduces the experiment performed by [START_REF] Plona | Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies[END_REF]. The models are made using the Discontinuous Galerking method (DGM) of [START_REF] Dupuy | Wave propagation in heterogeneous porous media formulated in the frequency-space domain using a discontinuous galerkin method[END_REF] and are checked against the reflectivity approach (SKB) used by De [START_REF] De Barros | Perturbations of the seismic reflectivity of a fluid-saturated depth-dependent poroelastic medium[END_REF]. The SKB synthetic data have been corrected from the 3D effects, using an infinite line of sources, in order to be directly comparable to the DGM solutions. A P-wave is generated with an explosive source (central frequency of 200 Hz) in an quasi-elastic layer (the Biot wave is entirely diffusive in this layer) and is transmitted and converted into S-wave and Biot wave at the interface with a porous layer. Receivers are set into the second layer and record the three waves. Fig. 3. Example of poroelastic and elastic numerical modelling to mimic differential data for CO 2 storage. Top) Resevoir models (left: Facies, right: Permeability) used in the modelling and modified from [START_REF] Manzocchi | Sensitivity of the impact of geological uncertainty on production from faulted and unfaulted shallow marine oil reservoirs -objectives and methods[END_REF]. Elastic properties are computed using the Gassmann relationships from the porous parameters. Middle) Left: Differential seismograms (model with CO 2 minus baseline) for elastic (red) and poroelastic data (black), and right: differences between the elastic and poroelastic differential data. Bottom) As the middle panels for a larger extension of CO 2 . Receivers and source are located on the free surface.

The explosive source has a 8Hz Ricker wavelet signature. (multiplied by a factor of 5). PP, PBiot and PS2 stand for the transmitted P and the converted Biot and S waves, respectively. PS1 stands for the conical wave associated with the direct P wave in the first layer. The explosive source is in a quasi-elastic layer, while receivers stand in a porous half space.

Two cases are studied: 1) a low frequency case, where the source frequency is smaller than the cut-off frequency ( f c = 6400 Hz). The Biot wave, in this case, is not propagative and cannot be seen; 2) a high frequency case, the source frequency is higher than the cut-off frequency ( f c = 0.64 Hz). This is obtained by decreasing the value of the fluid viscosity by 10000. The Biot wave becomes propagative, and can be observed, mainly in the fluid displacement data. In the seismograms of figure 4, the transmitted P-wave (PP), the conical wave associated with the P direct wave in the first layer (PS1), and the transmitted S-wave (PS2) are identified. These two S-waves (PS1 and PS2) will be uncoupled at further offsets. Finally, the very small differences between both modelling methods prove their accuracy.

A similar case is studied in a double porosity medium. Taking the same source/receiver layout with an explosive source in the first layer comprising a quasi-elastic sandstone and a line of receivers in the second layer comprising the double porosity medium described in the part 2.2 (sandstone with 3% of 1 cm spherical sand inclusions), we compute the seismograms and compare the double porosity results with the effective single porosity results. The seismograms of solid and relative fluid/solid displacements are given in figure 5. The influence of the double porosity homogeneization (via complex frequency dependent mechanic moduli) is clearly visible on the transmitted P-waves. Particularly, the waveforms are strongly distorted as the attenuation and dispersion are higher (see figure 2). As we are in the low frequency domain, the Biot waves are not visible in the seismograms, but they are responsible for the loss of seismic energy. As predicted by the theory, the converted S-waves are not impacted by the double porosity approach.

In these examples, we have demonstrated the importance of taking into account complex poroelastic theories in order to understand and reproduce real seismic signals whose waveforms are strongly impacted by the presence of fluid and medium heterogeneities.

Sensitivity analysis

In the next sections, we use the reflectivity algorithm SKB [START_REF] De Barros | Perturbations of the seismic reflectivity of a fluid-saturated depth-dependent poroelastic medium[END_REF] and focus on backscattered energy, i.e., we consider reflected seismic waves as in a seismic reflection experiment. We further assume that, whenever they exist, waves generated in the near surface (direct and head waves, surface and guided waves) are filtered out of the seismograms prior to the analysis. The assumption of plane-layered media is admittedly too simple to correctly describe the structural features of geological media, but it is nevertheless useful to explore the feasibility of an inversion process accounting for the rheology of porous media.

The sensitivity of the seismic waveforms to the model parameters is investigated for layered medium by computing the first-order derivatives of the seismic displacements with respect to the relevant poroelastic parameters. These operators, which are often referred to as the Fréchet derivatives, are expressed via semi-analytical formulae by using the Born approximation [START_REF] De Barros | Perturbations of the seismic reflectivity of a fluid-saturated depth-dependent poroelastic medium[END_REF]. They can be readily and efficiently evaluated numerically because they are only functions of the Green's functions of the unperturbed medium. In each layer, we consider the eight following quantities as model parameters: 1) the porosity φ, 2) the mineral bulk modulus K s , 3) the mineral density ρ s , 4) the mineral shear modulus G s , 5) the consolidation parameter c s , 6) the fluid bulk modulus K f , 7) the fluid density ρ f and 8) the permeability k 0 . This parameter set allows us to distinguish the parameters characterizing the solid phase from those describing the fluid phase. The fluid viscosity η is one of the input parameter but it is not considered in the inversion tests as its sensitivity is similar to the permeability.

Figure 6 (left) presents the sensitivity of the 8 parameters, for P and S waves. The sensitivity of the reflected wavefield varies drastically among the different parameters. We note that the reflected waves are especially sensitive to the mineral density ρ s , porosity φ, shear modulus G s and consolidation parameter cs. If we have some knowledge about the mineral properties (i.e., G s , K s and ρ s are fixed), the porosity φ and consolidation factor cs are the most sensitive parameters and therefore the key parameters to consider in an inversion procedure. On the other hand, the viscosity η and permeability k 0 have only a weak influence (10 4 times smaller than the porosity) on the wave amplitudes. The inversion for the parameters with such a low influence on the seismic waves will therefore be very delicate if other parameters are imperfectly known.

Figure 6 (right) shows the sensitivity of the P waves to the the fluid modulus K f and the mineral solid modulus K s . We note that the fluid modulus K f has a stronger influence than the solid modulus K s if the medium is poorly consolidated. The inverse is true for a consolidated medium. Similar patterns can be observed with the porosity: the higher the porosity, the stronger the influence of the fluid on the seismic waves. This means that it will be easier to determine fluid properties for an unconsolidated medium. For example, fluid substitution due to CO 2 injection leads to clear bright spot in Sleipner area (Norway, [START_REF] Arts | Monitoring of CO 2 injected at Sleipner using time-laps seismic data[END_REF], where the medium is poorly consolidated. The same set-up in stiff rocks does not produce such clear images, like in Weyburn field (Canada, White, 2009) . Left) Sensitivity of the reflected waves to the porous parameters for the P (downward triangles) and S (circles) waves, i.e., maximum of the energy reflected back by a perturbation of the porous parameters. Right) Sensitivity of the reflected waves to the fluid K f and the mineral bulk modulus K s as a function of the consolidation parameter cs. Note that higher consolidation parameter cs corresponds to softer materials.

To evaluate the coupling between parameters, we look at the Amplitude Versus Angle (AVA) curves in figure 7 for the PP (left) and SS (right) reflected waves due to a small and localized perturbation of a model parameter. We note that for some parameters, the model perturbations lead to similar modifications of the seismic response. For example, perturbations in densities and permeability show identical AVA responses. The same is true for the bulk moduli. This strong coupling between parameters will prevent simultaneous reconstruction of these parameters in an inversion process. [START_REF] Morency | Finite-frequency kernels for wave propagation in porous media based upon adjoint methods[END_REF] also investigated the sensitivity of the seismic waves in porous media. They determined finite-frequency kernels based upon adjoint methods and investigated different parameter sets, in order to find the set that leads to the minimal coupling between parameters. They concluded that decomposing the input parameters into seismic velocities is the most stable approach in an inversion code.

Full waveform inversion

Full waveform inversion has shown to be an efficient and accurate tool to study the subsurface in the acoustic and elastic wave theory [START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF]. Historically, most of the FWI methods [START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations[END_REF][START_REF] Tarantola | The seismic reflection inverse problem[END_REF] have been implemented under the acoustic approximation, for 2D model reconstruction (e.g. [START_REF] Gauthier | Two dimensionnal nonlinear inversion of seismic waveforms: numerical results[END_REF][START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF] or 3D structures (for instance, [START_REF] Ben-Hadj-Ali | Velocity model-building by 3-D frequency-domain, full waveform inversion of wide aperture seismic data[END_REF][START_REF] Sirgue | 3D frequency domain waveform inversion using time domain finite difference methods[END_REF]. Applications to real data are even more recent [START_REF] Hicks | Reflection waveform inversion using local descent methods: Estimating attenuation and velocity over a gas-sand deposit[END_REF][START_REF] Operto | Crustal-scale seismic imaging from multifold ocean bottom seismometer data by frequency-domain full-waveform tomography: application to the eastern nankai trough[END_REF][START_REF] Pratt | Seismic waveform inversion in the frequency domain, part II: Fault delineation in sediments using crosshole data[END_REF]. The elastic case is more challenging, as the coupling between P and S waves leads to ill-conditioned problems. Since the early works of [START_REF] Mora | Nonlinear two-dimensional elastic inversion of multioffset seismic data[END_REF] and [START_REF] Kormendi | Nonlinear waveform inversion of plane-wave seismograms in stratified elastic media[END_REF], the elastic problem has been addressed several times over the last years with methodological developments [START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF][START_REF] Choi | Two-dimensional waveform inversion of multi-component data in acoustic-elastic coupled media[END_REF][START_REF] Gélis | Two-dimensional elastic full waveform inversion using Born and Rytov formulations in the frequency domain[END_REF].

Using a poroelastic theory makes the problem even more difficult, especially because it adds much more unknowns. To the best of our knowledge, the first attempt to solve this Fig. 7. Energy of plane waves reflected from perturbations in ρ s , ρ f , k 0 , φ, K s , K f , G s , and cs, as a function of incidence angle. The eight upper panels and eight lower panels correspond to PP and SS reflections, respectively. The curves are normalized with respect to the maximum value indicated above each panel.

problem was made by De [START_REF] De Barros | Perturbations of the seismic reflectivity of a fluid-saturated depth-dependent poroelastic medium[END_REF] and De [START_REF] De Barros | Full waveform inversion of seismic waves reflected in a stratified porous medium[END_REF] for stratified media and by [START_REF] Morency | Finite-frequency kernels for wave propagation in porous media based upon adjoint methods[END_REF] and [START_REF] Morency | Acoustic, elastic and poroelastic simulations of co2 sequestration crosswell monitoring based on spectral-element and adjoint methods[END_REF] in 3-dimensional media. In the following sections, we will describe the main results obtained by De [START_REF] De Barros | Perturbations of the seismic reflectivity of a fluid-saturated depth-dependent poroelastic medium[END_REF] and De [START_REF] De Barros | Full waveform inversion of seismic waves reflected in a stratified porous medium[END_REF].

Inversion algorithm

Our method to determine the intrinsic properties of porous media is based on a full waveform iterative inversion procedure. It is carried out with a gradient technique to infer an optimum model which minimizes a misfit function. The latter is defined by a sample-to-sample comparison of the observed data d obs with a synthetic wavefield d = f (m) in the time-space domain, and by an equivalent term describing the deviations of the current model m with respect to an a priori model m 0 , i.e.,

S(m) = 1 2 ||d -d obs || D + ||m -m 0 || M , ( 6 )
where the L2-norms || . || D and || . || M are defined in terms of a data covariance matrix C D and an a priori model covariance matrix C M [START_REF] Tarantola | Inverse Problem Theory: Methods for data fitting and model parameter estimation[END_REF]. The model m contains the description of one or several parameters in layers whose thicknesses are defined by the peak content of the data [START_REF] Kormendi | Nonlinear waveform inversion of plane-wave seismograms in stratified elastic media[END_REF]. The model is updated using a quasi-Newton algorithm [START_REF] Tarantola | Inverse Problem Theory: Methods for data fitting and model parameter estimation[END_REF], which involves the Fréchet derivatives obtained earlier. As this problem is strongly non-linear, several iterations are necessary to converge toward an optimum model m, i.e, a model whose response d satisfactorily fits the observed data d obs .

Numerical results

In order to determine the accuracy of the inversion procedure for the different model parameters considered, we first invert for a single parameter, in this case the mineral density ρ s , and keep the others constant. The true model to reconstruct and the initial model used to initialize the iterative inversion procedure (which is also the a priori model) are displayed in figure 8. The other parameters are assumed to be perfectly known. Their vertical distributions consist of four 250 m thick homogeneous layers. Parameters φ, c s and k 0 decrease with depth while parameters ρ f , K s , K f and G s are kept strictly constant.

Vertical-component seismic data (labelled DATA, fig 9) are then computed from the true model for an array of 50 receivers spaced 20 metres apart at offsets ranging from 20 to 1000 metres from the source. The latter is a vertical point force whose signature is a perfectly known Ricker wavelet with a central frequency of 25 Hz. Source and receivers are located at the free surface. As mentioned previously, direct and surface waves are not included in our computations to avoid complications associated with their contributions. Figure 9 also shows the seismogram (labelled INIT) at the beginning of the inversion, i.e., the seismogram computed from the starting model. Figure 8 shows that the true model, which consists in 10 metre thick layers from the surface to 1000 metre depth, is very accurately reconstructed by inversion. As there are no major reflectors in the deeper part of the model, very little energy is reflected toward the surface, which leads to some minor reconstruction problems at depth. In figure 9, we note that the final synthetic seismograms (SYNT) almost perfectly fit the input data (DATA) as shown by the data residuals (RES) which are very small. The inversions carried out for the φ, ρ f , K s , K f , G s and c s parameters (not shown) exhibit the same level of accuracy. However, as predicted by De [START_REF] De Barros | Perturbations of the seismic reflectivity of a fluid-saturated depth-dependent poroelastic medium[END_REF] and [START_REF] Morency | Finite-frequency kernels for wave propagation in porous media based upon adjoint methods[END_REF] with two different approaches, the weak sensitivity of the reflected waves to the permeability does not allow us to reconstruct the variations of this parameter.

Being related to seismic wave attenuation and fluid flow, permeability appears as not only the most difficult parameter to estimate but also the one which would have the greatest benefits to the characterization of porous formations, notably in the oil industry (Pride et al., 2003). One possibility to estimate it is to measure the fluid motion, or, by reciprocity, to use fluid pressure sources.

As observed in the sensitivity study, parameters are strongly coupled. Multiparameter inversion is thus an ill-posed problem, which is, in most of the cases, not reliable, as errors on one parameter will map into the reconstruction of the other parameters. The use of analytical expressions for the sensitivity kernels allows an easy rearrangement of the parameter set, in order to invert for the most pertinent parameters. Using some a priori information, it is then possible to efficiently decrease the number of unknown parameters. For example, there is no reason to invert for both fluid parameters ρ f and K f , if we know that pores are filled by either gas or water. It is much more efficient to invert only for the saturation rate.

Differential inversion

To reduce the ambiguities of multiparameter inversion, a differential inversion has been considered and implemented [START_REF] De Barros | Full waveform inversion of seismic waves reflected in a stratified porous medium[END_REF]. Instead of dealing with the full complexity of the medium, we concentrate on small changes in the subsurface properties such as those occurring over time in underground fluid-filled reservoirs. This approach may be particularly useful for time-lapse studies to follow the extension of fluid plumes or to assess the fluid saturation as a function of time.

For example, the monitoring of CO 2 underground storage sites mainly aims at mapping the CO 2 extension. Time lapse studies performed over the Sleipner CO 2 injection site in the North Sea (see e.g. [START_REF] Arts | Monitoring of CO 2 injected at Sleipner using time-laps seismic data[END_REF] highlight the variations of fluid content as seen in the seismic data after imaging and inversion. In this fluid substitution issue, the parameter of interest is the carbon dioxide/saline water relative saturation. A differential inversion process will allow us to free ourselves from the unknown model parameters. This approach is valid for any type of fluid substitution monitoring problem, such as water-table variation, gas and oil extraction or hydrothermal activity.

The first step in this approach is to perform a base or reference survey to estimate the solid properties before the fluid substitution occurs. When performing a multiparameter inversion, the model properties are poorly reconstructed in general. However, the seismic data are reasonably well recovered. Thus, in spite of its defects, the reconstructed model respects the wave kinematics of the input data. In other words, the inverted model provides a description of the solid earth properties which can be used as a starting model for subsequent inversions. The latter would be used to estimate the fluid variations within the subsurface from a series of monitor surveys (second step).

To test this concept, we perform an inversion for two strongly coupled parameters, namely the porosity and the consolidation parameter. The resulting models are then used as starting models. We perturb the fluid properties of the true model to simulate a fluid variation over time. Two 30-metre thick layers located between 50 and 80 metres depth and between 110 and 140 metres depth are water depleted due to gas injection. The water saturation then varies between 60 and 80% in these two layers (figure 10). Our goal is to estimate the fluid properties by inverting the seismic data for the water saturation.

The model obtained is displayed in figure 10. We see that the location and extension of the gas-filled layers are correctly estimated. The magnitude of the water saturation curve, which defines the amount of gas as a function of depth, is somewhat underestimated in the top gas layer but is nevertheless reasonably well estimated. In the bottom gas layer, the inversion procedure only provides a qualitative estimate of the water saturation. These computations show that the differential inversion approach is capable of estimating, with reasonably good quality, the variations of fluid content in the subsurface without actually knowing all properties of the medium.

Conclusion

Using a poroelastic theory is much more complex than an elastic or a visco-elastic theory. However, poroelastic theories are an attempt to quantitatively describe the attenuation processes from the physical properties of the geological material. Furthermore, at seismic frequencies, attenuation is dominated by the mesoscopic scale mechanism, involving fluid flow at the boundaries of any heterogeneities. Poroelastic theories intrinsically take into account this loss of energy, while the equivalent visco-elastic approach neglects it. As the near surface media are stronly heterogeneous, with strong lateral and vertical contrasts, and different fluids involved, one has to deal with full poroelastic theories to accurately consider attenuation and fluid-solid motions.

The sensitivity to the different parameters varies hugely among parameters, and parameters are strongly coupled. Using a poroelastic theory to reconstruct model properties is in its nature an ill-conditioned problem. It shows however very promising possibility for differential inversion, and for certain issues where the problems can be reduced to the determination of only few parameters.

Poroelastic theories are, of course, not perfect yet, as they fail to give an universal law to explain seismic wave attenuation and propagation. They are however the direction to go if one wants to use full waveform inversion to make quantitative imagery of the rock physics and subsurface fluids. In particular, the permeability is a key parameter for exploration; the reconstruction of such a parameter from seismic waves will necessitate the use of complex poroelastic theories. Development in this direction still has to be continued. In particular, for imagery problems, data have to be improved to get around the problems of coupling between parameters. This can be done by using the information carried by fluid motions, which can give new insights into the permeability and the fluid properties, or by exploring in deeper details the coupling between seismic and electromagnetic waves.

  Fig.1. Example of data for which a poroelastic-based interpretation may be useful. Data are recorded on a seashore in the South of France. Source is a hammer shot, and the 24 receivers are equally spaced between 10 and 55 m from the source. The medium is very soft and water saturated, with a direct P wave velocity of c. 1600m/s and a S wave velocity lower than 200 m/s.

  Fig.2. P-wave velocities (top panel) and attenuations (bottom panel) with respect to the frequency (between 1 and 10 8 Hz). Blue lines are for the double porosity theory for inclusion sizes a = 1 cm (continuous lines), a = 5 cm (dotted-dashed lines) and a = 10 cm (dashed lines). Poroelastic responses of the individual phases are the black lines. The results for the average single porosity theory, computed using weighted average of each phase properties, are given by the red lines.

Fig. 4 .

 4 Fig. 4. Flat interface case: Left) Seismograms of (top) vertical solid u z , and (bottom) relative fluid/solid displacement components w z in the low-frequency regime; Right) As Left, but in the high-frequency domain. The SKB solution is indicated by a continuous line and the DGM by crosses; dashed-dotted lines indicate the differences between the two solutions (multiplied by a factor of 5). PP, PBiot and PS2 stand for the transmitted P and the converted Biot and S waves, respectively. PS1 stands for the conical wave associated with the direct P wave in the first layer. The explosive source is in a quasi-elastic layer, while receivers stand in a porous half space.

Fig. 5 .

 5 Fig.5. Double porosity case: seismograms of top) solid displacement u x (left) and u z (right) and bottom) relative fluid/solid displacements w x (left) and w z (rigth). The setup is the same as for figure4, the explosive source (Ricker with a 200Hz peak frequency) is in a quasi-elastic layer, while receivers stand in a double porosity half space. The double porosity solution is indicated by a continuous black line and the effective single porosity by a dashed blue line. PP and PS stand for the converted P and S waves, respectively.

  Fig.6. Left) Sensitivity of the reflected waves to the porous parameters for the P (downward triangles) and S (circles) waves, i.e., maximum of the energy reflected back by a perturbation of the porous parameters. Right) Sensitivity of the reflected waves to the fluid K f and the mineral bulk modulus K s as a function of the consolidation parameter cs. Note that higher consolidation parameter cs corresponds to softer materials.

Fig. 8 .Fig. 9 .Fig. 10 .

 8910 Fig. 8. Models corresponding to the inversion for the mineral density ρ s : initial model, which is also the a priori model (dashed line), true model (thick grey line), and reconstructed model (black line). The corresponding seismograms are shown in figure 9.

Using a Poroelastic Theory to Reconstruct Subsurface Properties: Numerical Investigation

Seismic Waves, Research and Analysis

Seismic Waves, Research and AnalysisUsing a Poroelastic Theory to Reconstruct Subsurface Properties: Numerical Investigation

Acknowledgements

L. De Barros and G. S. O'Brien were partly funded by the department of Communications, Energy and Natural Resources (Ireland) under the National Geosciences programme 2007-2013. B. Dupuy, S. Garambois and J. Virieux are supported by the National Research Agency (ANR) "Captage et Stockage de CO 2 " program (ANR-07-PCO2-002). The numerical computations were performed by using the computational facilities from SFI/HEA ICHEC (Ireland), the HPC national computer centers (France, CINES and IDRIS under the allocation 2010-046091 GENCI) and the HPC center of the Grenoble observatory (OSUG).


Cyril Buttay
email: cyril.buttay@insa-lyon.fr

Amandine Masson

Jianfeng Li

Mark C Johnson

Mihai Lazar

Christophe Raynaud

Hervé Morel

Die Attach of Power Devices Using Silver Sintering -Bonding Process Optimization and Characterization

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

Introduction

Soldering is the most common die-attach technique for power electronic devices. Indeed, such components require a bond that is both electrically and thermally conductive. Classical solder alloys include lead-free SACs ( Sn, Ag, Cu, with a melting temperature slightly higher than 200°C), HMP (High Melting Point, with a high lead content, around 300°C), and gold-based alloys (Au-Sn, Au-Ge, AuSi, up to 362°C).

High temperature electronics obviously require die attaches that have a high melting point. Furthermore, as a material approaches its melting temperature, its mechanical performances decrease dramatically. Therefore, a safety margin should be kept between the operating temperature and the melting point of the die attach material.

The ratio between the operating and melting temperatures (expressed in kelvin) of a material is called the homologous temperature. For electronics applications, it is commonly as high as 0.8 or more. However, even with such a high value, it can be seen that as the operating temperatures increase, soldering alloys with much higher melting points are required: for T operating = 200°C (473 K), a melting point of 473 0.8 = 591 K (318°C) is necessary. For T operating = 300°C, the solder must remain solid up to 443°C. This requires a high temperature bonding process, resulting in an increased level of residual stresses.

A survey of high-temperature die attach solutions has been recently published [START_REF] Manikam | Die Attach Materials for High Temperature Applications: A Review[END_REF]. Among the solutions presented, some are not based on soldering, i.e the bonding process can be performed at a temperature lower than the melting point of the resulting joint. They are attractive, as they present the advantage of a low-temperature process (low level of residual stress) and of the high melting temperature (wide safety margin between operating temperature and melting point, or low homologous temperature value).

The first solution is the Transient Liquid Phase Bonding (TLPB), in which a thin layer of a low-melting point material A (tin, indium. . . ) acts as a liquid phase between two thicker solid layers of material B. As material A diffuses into B, isothermal solidification occurs. Eventually, one ends up with a single phase, with a melting temperature higher than that of material A. Results have been presented using a thin layer of tin between copper layers [START_REF] Guth | New assembly and inter-Figure 10: FIB ion channelling contrast image of a silver joint between a Si die and an DBC substrate (FIB micro-section) connects beyond sintering methods[END_REF] or between silver layers [START_REF] Mcpherson | Packaging of High Temperature 50kW SiC Motor Drive Module for Hybrid-Electric Vehicles[END_REF]. In some cases, indium was used instead of tin [START_REF] Quintero | High Temperature Die Attach by Transient Liquid Phase Sintering[END_REF]. Other solutions exist, based on gold-tin, bismuth, . . . Another approach is based on powder sintering. This technique has been studied for the last 20 years [START_REF] Schwarzbauer | Method of securing electronic components to a substrate[END_REF]. A metal powder is mixed with organic materials (thinner, binder, dispersant) to form a paste. Silver is usually chosen as the metal, because of its excellent thermal and electrical performances [START_REF] Schwarzbauer | Novel large area joining technique for improved power device performance[END_REF], its reasonable price and its low temperature sintering capability. Furthermore, sil-ver can be sintered in the presence of air. The silver paste is printed or dispensed on a substrate, exactly like a solder paste, and the die is placed on top of the paste. Pressure (several tens of MPa) and temperature (around 240°C) are then applied for the sintering to occur. Such die attach solution is being used for industrial production [START_REF] Göbl | Low temperature sinter technology Die attachment for automotive power electronic applications[END_REF].

The pressure applied during sintering requires specific tooling to avoid breaking the semiconductor dies. To avoid or reduce the need for pressure, it is possible to use silver particles with a smaller size. Indeed, reducing the particles size results in a higher specific surface (more particle surface per powder volume), and makes them more prone to sinter. This can be used either to reduce the sintering temperature (room temperature sintering was demonstrated in [START_REF] Wakuda | Ag Nanoparticle Paste Synthesis for Room Temperature Bonding[END_REF]) or to reduce the sintering pressure [START_REF] Bai | Low-Temperature Sintering of Nanoscale Silver Paste for Semiconductor Device Interconnection[END_REF]. Zero-pressure sintering has been demonstrated, but it was also found that applying a reasonable amount of pressure (5 MPa) improves the bond quality, especially for the largest dies (1 cm 2 ) [START_REF] Lei | Low-Temperature Sintering of Nanoscale Silver Paste for Attaching Large-Area (>100mm 2 ) Chips[END_REF] In the present article, we will focus on this nanoparticle silver paste. The objective is to study the bonding parameters, in order to achieve a well controlled (i.e. reproducible) die attach, and to describe the process details.

2 Die attach procedure

Silver Paste

The paste selected is the Nano Ag paste from NBE Tech [START_REF] Bai | Low-Temperature Sintering of Nanoscale Silver Paste for Semiconductor Device Interconnection[END_REF], containing additives (binder, dispersant, thinner) and nano Ag particles (8-12 nm size). Several articles have focused on the bonding process using this paste. The main parameters being studied were the temperature profile [START_REF] Wang | Low-Temperature Sintering with Nano-Silver Paste in Die-Attached Interconnection[END_REF] (number and temperatures of the steps, ramping time. . . ) and the die size [START_REF] Lei | Low-Temperature Sintering of Nanoscale Silver Paste for Attaching Large-Area (>100mm 2 ) Chips[END_REF]. Therefore, the temperature profile used in this study (Fig 1) could be directly selected from the literature.

Substrate and die

The main substrate used for this study is Direct Bonded Copper (DBC) from Curamik (635 µm of Al 2 O 3 with 300 µm of copper clad on both sides). For some tests, one of the copper surface was polished in order to reduce the surface roughness (see Fig 2). In both cases, the copper surface was cleaned immediately before the assembly: 5 minutes dip in H 2 SO 4 (96%) to remove oxidation, followed by rinsing in deionized water, and 3 minutes in ethanol (US bath). The samples were then dried using a nitrogen blow nozzle.

For two series, other substrates were used: one series was made using a Curamik DBC with a nickel-gold finish and another with a Denka substrate (635 µm Si 3 N 4 with 300 µm copper), also with a nickel-gold finish. Although the gold layer provides protection against oxidation, the cleaning process was kept identical (H 2 SO 4 , rinsing, ethanol, drying).

The SiC dies are mechanical samples ("dummy chips") provided by SiCED. They measure 2.7 × 2.7 mm, with a silver backside finish. Prior to assembly, they were cleaned as follows: 3 minutes in ethanol (US bath), followed by drying using a nitrogen blow nozzle.

Application of the Paste

The silver paste can was first opened in april 2009, and has since been stored at room temperature. Although it Figure 3: The manual screen-printing system, with its stainless-steel squeegee. The screen size is roughly equivalent to the A4 paper format is supposed to have a pot life of one year only, no change in consistency or in sintering behaviour was noticed. It was therefore decided to pursue the tests with the same paste, even once it has passed its expiry date.

The paste was applied by screen printing, to ensure a uniform deposit thickness. Two stencils were used: 50 µm and 100 µm. While paste manufacturers advise to use the thinner, it is much more expensive (around 500 C) than the thicker (around 150 C). This is because thin screens are obtained by electroforming (similar to electroplating), while screens of 100 µm and thicker are laser cut in a metal foil.

The screens were attached to a simple frame, with a manual squeegee (fig. 3). The paste has a thick consistency, making it relatively easy to print, with a good shape retention. However, if it is simple to print, the paste is very difficult to clean. Extra care should be taken not to stain any surface with the paste. Nickel or stainless steel screens can be cleaned without too much trouble (using isopropyl alcohol or ethanol), this is not the case with more porous surfaces such as anodized aluminium. In many cases, cleaning attempts simply results in the spreading of the paste! Furthermore, as the effect of silver nanoparticles on living cells is still investigated [START_REF] Lubick | Nanosilver toxicity: ions, nanoparticles or both?[END_REF], it is important to wear gloves when handling stained equipments.

Sintering

As described in section 1, the main objective of using nano-particles of silver is to reduce the pressure required to less than a few mega-pascals. Most presses are designed to apply forces upwards of a kilo-newton, and the silicon carbide dies tend to be small (in our case 2.7 × 2.7 mm). This results in pressures exceeding 100 MPa. Therefore, a specific press was developed for our application. A photo can be seen in fig. 4 and the mechanical structure is described in fig. 5.

This press can apply forces ranging from 5 N to approximately 1 kN. This latter value is require to attach several dies at the same time. Both heating platens can be controlled synchronously or independently, and the pressure can be applied or removed at any time during the bonding cycle. Both platens can achieve 50°C/min ramping time and are equipped with 750 W heating power (cartridge heaters for the top platen, lamps for the bottom platen). A force sensor has recently been installed (not visible in fig. 4, but represented in fig. 5) to measure the actual force applied by the press. This force can be adjusted by changing the air pressure in the pneumatic cylinder.

The parallelism between both platens is ensured by a moving frame which slides (using ball bearings) along four shafts. The top platen is secured to this moving frame via 4 springs in order to accommodate a possible tilt of the die. To bond several dies simultaneously, or if the die is fragile (large and thin) a soft material (viton™ rubber or PTFE) is placed between the dies and the top platen. On the contrary, to bond a single, small die, no interface material is used, the die is in direct contact with the top aluminium platen.

After paste printing, the substrate is put in the press (fig. 4). The die is then carefully placed on top of the paste deposit. No force is applied on the die at this stage. The press is then turned on.

For most of the assemblies described in this paper, the pneumatic cylinder (see fig. 5) is actuated from the beginning of the cycle in fig. 1. This ensures a good contact between the die and the fresh paste.

For some assemblies, however, a two-step printing of the paste was used. This method is described in [START_REF] Lei | Low-Temperature Sintering of Nanoscale Silver Paste for Attaching Large-Area (>100mm 2 ) Chips[END_REF]: the paste is printed on the substrate, and dried at 150°C without any pressure to remove the organic materials. Then a fresh lump of paste is printed, using the same screen, and the die is placed on top of the deposit. The rest of the process is identical to the one-step technique (temperature profile of fig. 1). This technique is supposed to give better results with large dies.

Effects of the process parameters

A first series of experiments is listed in table 1. For each series, at least 6 die attaches were performed: 1 was kept for micro-section and microscope inspection, while the remaining assemblies were used for mechanical (destructive) testing. The objective was to find the parameters that would not only result in a strong bond, but more importantly that would result in a reproducible bond. That is, to find the parameters for which all dies attaches would break for the same force applied. Due to the shape of the assemblies, the most practical (and commonly used) bond strength assessment is the so-called "shear test" in which a tool applies a horizontal force on one side of the die, the substrate being kept secured. This results in a shear stress in the die attach. The shear strength is the maximum force recorded before the failure of the die attach divided by the die surface area. The test equipment used is visible in fig. 6.

The results are visible in figure 7. It can be seen that for most series, the results are spread over a wide range.

Only series E show an acceptable level of reproducibility.

Two mechanisms can explain the poor reproducibility of all series except series E:

• The roughness of the substrate: as can be seen in figure 8 (as well as in figure 2), the roughness of a "raw" DBC is comparable with the thickness of the silver joint. This means that on some places, the die can be in direct contact with the DBC, while on other, the silver layer is not thick enough to prevent voiding. On a polished DBC, on the opposite, the silver joint is of regular thickness. This could explain the spreading of series A', B, C, D, F, J

• The oxidation of the copper layer: for the two-step deposits (series N and O), the copper under the silver deposit has more time to oxidize, resulting in a weaker silver-to-copper bond.

As a result, the most satisfying parameter set was found to be: 50 µm of silver paste on a polished substrate, with 6 MPa pressure applied on the die. Note that the silver layer thickness after sintering was found to be under 10 µm, which is much less than expected (around half the thickness of the fresh paste, or 25 µm). This is because a large part of the paste is "squeezed" out when applying the pressure. In series F, a thicker screen was used to deposit 100 µm of paste, but no improvement was observed. The pressure was reduced for series J, but that resulted in an even less predictable result. In series N and O, the two-stage deposit (also aimed at producing a thicker silver layer) did not show any improvement over series A'.Finally, series B, in which a manual scrubbing was performed on the dies to improve their contact with the silver paste, show a very poor adhesion. This can be explained by the scrubbing motion, which tends to further reduce the silver layer thickness. As the failure mechanism for all the assemblies was a fracture along the silver-substrate interface, improving the die-silver interface did not bring any improvement.

For all assemblies performed on bare copper substrates, it was found that the exposed copper surfaces were deeply oxidized. This means that a gold-finish is desirable (silver sintering was found not to adhere directly on nickel). The substrates used for series D only have a gold-flash finish (50 nm) on a nickel plate, but this thin gold layer was found to be sufficient to record the best shear strength (52.3 MPa).

Other experimental considerations 4.1 Simultaneous bonds

With the chosen temperature profile, which is fairly slow compared to soldering (almost 3 hours including cooling-down), it is important to bond several dies at once. As the silver layer is very thin, it does not allow to accommodate any difference in height between dies, or any irregularity in substrate thickness. Therefore, an interface material is necessary between the dies and the top platen or the press. Two materials were assessed: PTFE and Viton. Both are 1.5 mm thick. PTFE is more rigid, while Viton has a rubber-like consistency.

With both materials, a satisfying bond could be achieved. However, the dies left a deep imprint in both interface materials. Also, the silver paste "squeezed out" from under the die stuck on the interface material, meaning that it could not be reused without polluting the top of the following dies with silver. This is an issue, as both materials are relatively expensive.

Large dies

9 × 9 mm and 12.7 × 12.7 mm silicon dies were used to evaluate the practical aspects of bonding large dies. As with the smaller SiC dies, the first attempts were conducted with a fairly stiff interface material between the dies and the top platen (a graphite foil).

As it can be seen in figure 9, this resulted in the bowing of the die: the silver layer is found to be much thinner at the edges of the die than in the middle. This is because the top of the die is not flat: the edge passivation (several tens of microns thick) protrudes from the die and makes contact with the platen, while the top metallization is in a recess and does not touch the platen. Therefore, even for a single die, a soft interface material is required for large dies to approach quasi-isostatic conditions.

The image of one of the silver layers is visible in figure 10. This image shows that what appears to be a solid layer of silver with the optical microscope actually consists of 200 to 1000 nm Ag grains with 150 to 400 nm pores/voids.

Conclusion

The implementation of a silver-sintering based solution for die attach has been described. Although this technology does not require very complex equipments (the process takes place in air, with moderate temperature and pressure), some precautions must be taken to ensure good results.

In particular, it was shown that the roughness of the substrate has a very significant effect on the reproducibility of the die attach. Other effects, such as substrate finish or the thickness of the paste deposit, are not as obvious. In particular, further experiments are needed to reduce the pressure required and achieve thicker silver joints.

Figure 1 :Figure 2 :

 12 Figure 1: The temperature profile used for the single-step sintering

Figure 4 :

 4 Figure 4: This sintering press (specific developpement) is equiped with two heating platens (heating rate exceeding 50°C per minute) and can apply forces ranging from a few newtons up to several hundred newtons.

Figure 5 :

 5 Figure 5: mechanical structure of the press. The mobile frame is attached to two counterweights, to achieve a low minimal force (<5 N).

Figure 6 :

 6 Figure 6: The shear tester used (specific development), with a 20 kgf (200 N) load cell mounted

Figure 7 :

 7 Figure 7: Shear strength measured for the series described in table 1. The line in each box corresponds to the median value (2 nd quartile), and the lower and upper limits of the box represent the first and third quartile respectively. The "whiskers" correspond to the extreme values. In other words, 50 % of the samples are in the box, 25 % in the lower whisker and 25 % in the top whisker.

Figure 8 :

 8 Figure 8: cross section of a sample of series A' (raw DBC, (a)) and of series E (polished DBC, (b)). It can be seen that the roughness of the raw DBC is comparable to the silver joint thickness.

Figure 9 :

 9 Figure 9: cross section of the die attach of a 9 mm Si die at the edge (a) and in the middle (b).

Table 1 :

 1 Test matrix for the SiC dies

	Series	thickness	substrate	bonding pressure remarks
	A'	50 µm	Raw DBC	6 MPa
	B	50 µm	Raw DBC	6 MPa	manual scrubbing to get good die-paste contact
	C	50 µm	Au-finished DBC	6 MPa
	D	50 µm	Au-finished Si 3 N 4	6 MPa
	E	50 µm	polished DBC	6 MPa
	F	100 µm	Raw DBC	6 MPa
	J	50 µm	Raw DBC	0.7 MPa
	N	50+50 µm	polished DBC	6 MPa	two-step screen printing
	O	50+50 µm	Raw DBC	6 MPa	two-step screen printing

Acknowledgments

The authors would like to thank the "Fondation pour la Recherche en Aéronautique et Espace" (FRAE) its financial support of this work, under grant name "EPAHT"


Dainis Zeps
email: dainize@mii.lu.lv

On building 4-critical plane and projective plane multiwheels from odd wheels

Keywords: graph coloring, chromatic critical graphs, wheels, planar graphs, projective planar graphs, Grötzsch graph, Mycielski's construction

.

Introduction

We are using terminology from [START_REF] Kratochvíl | About minor closed classes and the generalization of the notion of free-planar graphs[END_REF][START_REF] Bojan | Graphs on Surfaces[END_REF][START_REF] Zeps | Free Minor Closed Classes and the Kuratowski Theorem[END_REF].

... ... ... ... . We consider graph in fig. 1 left, that in [START_REF] Zeps | 4-critical wheel graphs of higher order[END_REF] is denoted G 3 , but in this article gets several denotations due to its particular features, w 111 or w 13 , g 111 or g 13 , see lower. In this introduction we consider some simple features of this graph w 13 called base graph in frames of this article, but further in the article we generalize these features to two, plane and projective plane, unbounded classes of graphs.

We start from a simple observation that contraction of one of edges (incident to degree 3 vertices) in the graph (w 13 ) turns it into octahedron graph (O) minus an edge (O -), see fig. 1. Now it is obvious that the graph O -is a minor of w 13 , but the graph O is not.

Further, we may express this fact using terminology from [START_REF] Zeps | 4-critical wheel graphs of higher order[END_REF], i.e., < O -, O > is minor bracket for the graph w 13 . We remind that < h 1 , h 2 > is minor bracket for G if h 1 ≺ h 2 , h 1 ≺ G and h 2 ≺ G, [START_REF] Zeps | 4-critical wheel graphs of higher order[END_REF]. The fact expressed as predicate < h 1 , h 2 : G >.

One more observation that the base graph w 13 turns into the wheel graph W 3 after contracting three edges incident to three and four degree vertices.

Further, graph w 13 is 4-critical [START_REF] Zeps | 4-critical wheel graphs of higher order[END_REF], and this fact might be verified directly for such a small graph. But for further discussion we need to examine this graph more closely.

The base graph w 13 can be considered as edge sum modulo two of three graphs W 3 in the way that each pair overlap just in one edge, and all three wheels have one vertex in common, see fig. 2. Let us assume that overlap either two rim edges, or two spike edges, or spike edge and rim edge, all three ways giving the same graph w 13 (see fig. 2): trivial (automorphism) fact, but not so for further, see below. Further we are going to use this consideration of w 13 as sum of three wheels, and for that reason we use this multi-index denotation for it, i.e., w 111 (or w 13 ) with three indices (three ones), not single index, and second index as factor of equal indices, see lower. Further we are going to use by edge summation modulo two denotations with simple arithmetic operators, i.e., multiplication and sum operators, e.g. w 13 = 3W 3 . Of course, corresponding graph operations behind are indeterministic, i.e., depending on how we configure graphs one against other by edge summation modulo two, i.e., which elements of corresponding graphs we allow to overlap. Besides, w 13 may be expressed in the way O + W 3 with rim edges annihilating with a triangle of octahedron graph, thus having one more equation w 13 = O + W 3 becoming equality under specified conditions.

The graph w 13 can be colored in four distinct ways, see fig. 3. We distinct them in way lonely color, D, is applied: the lonely D may color central hub, section hub, rim and two vertices of rim, thus giving four ways. We remind that only in chromatic critical graph every vertex may receive lonely color. Thereby, each vertex for w 13 may be colored with lonely color, thus proving that graph is 4-critical.

Further, the base graph w 13 may be embedded on projective plane, quadrangulating it, see a) in fig. 8 such on projective plane. It is remarkable that dual graph to octahedron is cube graph, but cube graph with one corner cut off becomes just the graph w 13 , that is selfdual (see [START_REF] Zeps | 4-critical wheel graphs of higher order[END_REF]). The aim of this article is to show that these facts concerning this single and very simple graph w 13 may be extended for unbounded classes of 4-critical graphs, both planar and projective planar. These classes we are going to build from arbitrary odd wheel graphs.

The theorem about octahedron minor bracket is mentioned in [START_REF] Zeps | Application of the Free Minor Closed Classes in the Context of the Four Color Theorem[END_REF], but here we give proof of this fact.

Defining 4-critical plane multiwheels

Let us assume that the wheel W k is built from k simple sections where as simple section we take triangle C 3 with one vertex common from each triangle for the wheel's hub, and opposite edge of triangle as for forming wheel's rim. In other words, wheel W k is the sum of its k sections (triangles) into a simple graph (with double-edges turning into simple edges). By the way, if we applied this same summation of the edges modulo two then the sides of triangles that touch each other would annihilate, and remnant/resulting graph would be simple cycle (rim of the wheel) and isolated vertex (hub of the wheel). As an obvious observation, let us notice that edge sum modulo two of k > 3 triangles gives C k and isolated vertex only in case triangles annihilate two edges, and one vertex becomes common for all triangles.

If we take wheel W k to be odd (k is odd) then we get 4-critical graph. Odd wheels are simplest unbounded class of vertex-3-connected 4-critical graphs. (Let us remember that odd cycles is the only graphs that are 3-critical.) We are going to generalize this class to 4-critical vertex-3-connected plane and projective plane multiwheels.

Let us do similar as higher summation of edge sets of some wheels. Let us replace each section in a simple odd wheel with another arbitrary odd wheel in a way that edge sets of sections are summed modulo two. This new aggregation of wheels M may be expressed as W ki where summation is modulo two over index i numbering sections that were replaced by wheels of order k i in each case. Evidently the summation itself is indeterministic because result of it depends on how wheels overlap each other by summation. Now we ask: under which conditions the resulting graph is 4-critical. It turns out that the answer directly is connected with the number of edges wheels intersect by summation. Theorem 5 below says that this number must be equal to two.

But first we are to prove some lemmas about intersecting wheels without edge losses by summation modulo two. For example, W 5 + W 5 may be formed with all five rim vertices of both wheels common and forming subgraph K 5 , but resulting 6-chromatic graph is in no way chromatic critical.

Lemma 1 Let k odd wheels by summation possibly intersect in vertices but not in edges. The sum of these wheels can't give 4-critical graph except in case k = 1.

Proof. Let by summation of wheels k > 1 wheels have become vertex connected components, that are wheels all the same, in the resulting graph. It is obvious that elimination of edge or vertex, or edge contraction in one wheel can't affect coloring of others in four colors in the resulting graph. Thus, the graph can't be chromatic critical.

Let us configure two odd wheels so that they have common two adjacent vertices, and sum their edges modulo two. It is easy to see that the resulting graph is 3-chromatic. It suffices to notice that losing of an edge in both odd wheels allow to color them in 3-colors so that lost edge's ends receive the same color. Further, we may easily apply the use of this fact to unclosed sequence of wheels W q1 ,...,W qi where two proximal wheels overlap in two adjacent vertices, but next two possibly only in one. Let us formulate it as a lemma.

Lemma 2 Let summation of edges modulo two is applied to unclosed sequence of wheels. The resulting graph is 3-chromatic.

We need one more crucial feature of 4-critical graphs. Let graph H is 4-critical and let H ′ = H ⊙ w be graph H with vertex w split into two new vertices and edges incident to w be connected either to one or other vertex. We ask whether graph H ′ can remain to be 4-critical. Of course, it is expectable that H ′ becomes 3-chromatic.

Lemma 3 Let graph H be 4-critical and w ∈ V (H). Then graph H ⊙w is always 3-chromatic.

We are not going to proving this fact here, but delegate to article [START_REF] Zeps | Can graph remain k-critical after vertex split?[END_REF]. We call the feature of 4-critical graph to loose his state of 4-critical graph under vertex split preservation of 4-criticality under vertex split. Thus, k-criticality is the feature that is preserved under vertex and edge elimination and edge contraction and expectably under vertex split too. It is easy to see that 3-criticality is preserved under vertex split and we ask this same for arbitrary k-critical graphs.

We are going to use lemma 3 in the following way. By summation of edge sets of wheels modulo two in order to build new 4-critical graphs we may ignore cases where wheels intersect only in vertices without incident edges, knowing that this can't lead to new 4-critical graphs. Suppose we received 4-critical graph in this way. Then splitting all vertices that were merged by summation backwards we should receive 3-chromatic graph but it might not be true. Let us formulate this fact as lemma.

Lemma 4 Let by summation of edge sets of wheels some wheels intersect in vertices without incident edges. Then resulting graph can't be 4-critical. Now we may go over to the main theorem of this chapter.

Theorem 5 Let 2k + 1 (k > 0) arbitrary odd wheels be summed in a way that edges of wheels are summed modulo two and all wheels have one overlapping vertex. The resulting graph M is 4-critical if and only if each wheel by summation modulo two looses just two of its edges and resulting graph is planar.

Proof. Let us first observe that two wheels may overlap in one or three edges but not in two, thus, to get two annihilating edges a wheel should overlap with two other wheels with one overlapping edge in each.

Let us first assume that wheel graphs are summed observing two edge loss condition. In this case four configurations of new section are possible, see fig. 4. We denote graphs achieved in this way by w k1-k2-...kq where q is number of sections, where in each section there are 2k i + 1 edges, and call them multiwheels. In this type of denotation we as if ignore three/four ways of section's configuration, but one could easily elaborate denotation with taking these different types of sections into account, see below.

In figure 2 we see simplest case where summed are three wheels W 3 giving multiwheel w 1-1-1 . In place of multiindex k 1 -k 2 -...k q with hyphens we equally use denotation without hyphens in case no confusion might arise, e.g. w 111 in place of w 4. Four types of sections for planar multiwheels possible. If third type we take in two oriented ways, leftwards and rightwards, then we get one type of section more. When so numbered first two and two other sections are mutually dual. It is convenient to characterize type of section by type (rim or spike edge) of lost edges in wheel. Then sections are 1) rim-rim-section (rr-section), 2) spike-spike-section (ss-section), 3) spike-rim-section (sr-section) and 4) rim-spike-section (rs-section).

It is easy to see that w 111 is 4-critical. Indeed, reserving central hub vertex as eventual lonely color vertex (receiving color D, see fig. 3) other vertices get forced colors. Further, both hub edge, spike edge and rim edge contractions/eliminations lead to 3-chromatic graph.

Let us assume that the resulting graph is planar. In that case it is convenient to characterize type of section by type (rim or spike edge) of lost edges in wheel. Then sections are 1) rim-rim-section (rr-section), 2) spike-spike-section (ss-section), 3) spike-rim-section (sr-section) and 4) rim-spike-section (rs-section), see fig. 4. Section arisen from w 1 we call simple section, which is of arbitrary type due to automorphisms, i.e., rim edges are spike edges too, and reversely. Now let us consider first type of section, rr-section, fig. 4. At least one vertex on rim may receive third color and then corresponding section hub receives forth color. Removing rim edge makes possible to color rim with two colors, but removing section's spike edge allows now both previous adjacent vertices color with one color, thus avoiding fourth color.

Let us consider second type section. Now the same applies for the local rim edges of the section. At least one vertex of the local rim should receive third color, and corresponding local rim edge or spike edge elimination may avoid use of fourth color.

Let us consider third (and fourth) type of section. Now outer hinges and vertex adjacent to central hub should receive different colors, but removal at least one edge from section violates this condition and allows to color hinge vertices of inner rim with the same color.

Thus, we have proved that multiwheel is 4-critical.

Let us prove theorem in the other direction.

According lemmas 1,2,4 we are to consider only those sums of wheels where intersections of vertices without incident edges are absent and unclosed sequences of wheels are absent. Even more, if some closed sequences are present, but some wheels as unclosed ends are present, these cases are not to be considered be-cause can't give 4-critical graphs. The only cases are these where only closed sequences of wheels are present. Further, only one closed sequence as cycle is to be considered for further.

Further, let us consider case of non-planar resulting graph, see fig. 5. In that case we have engaged in cycle of wheels some with two spikes that are not sequencing one to other. In these cases such wheel may be as if in two ways taken in cycle along one or other orientation in this wheel between two spikes. It is easy to see that this leads to fact that such wheel cant't give resulting graph as 4-critical, because odd wheel may be divided only in half of odd wheel and half of even wheel. Taking into cyclic path even part of wheel would spoil oddtimes-odd structure necessary for 4-critical multiwheel. (See fig. 5, where rights minimal case of non-planar resulting graph is shown to be 3-chromatic.) Thus, we can't afford non-planar spike pairs in edge summation modulo two. We have come to conclusion that closed sequence should be planar. Fig. 5. Illustration to proof. Minimal possible case of nonplanar section for as if eventual nonplanar multiwheel w115. Left, we see non-planar section, where annihilated spikes (dotted) are not sequencing ones, and section falls into as if two subwheels, one even subwheel with subrim 1-2-3 and one odd subwheel with subrim 3-4-5-1. Right, we color this nonplanar multiwheel into three colors, i.e., it isn't 4-critical, even not 4-chromatic.

Thus, under assumption that all summing wheels have in common one vertex each in each wheel with two incident edges that annihilate by summation modulo two we have proved what was necessary. We have come to 4-critical graph only by specified conditions.

We have proved the theorem.

Denotations for plane multiwheels. Some characteristics

Let us introduce some notational conventions for wheels and multiwheels. Along with traditional denotation for wheels with capital letter W with index k, i.e. W k , for odd wheels of order k = 2q + 1 we use denotation w q . For plane multiwheels we use letter w with odd (k ≥ 3) indices, w q1,...,q k where w qi was i-th wheel in summation. Let us introduce quantity Q = q i and sometimes use denotation w Q for this multiwheel.

If odd wheel has order k = 2q + 1, it has n = 2q + 2 vertices and m = 2n -2 = 4q + 2 edges. Similar expressions hold for multiwheels, where in place of q stands Q. Indeed, multiwheel has 2Q + 1 vertices and 4Q or 2n -2 edges. To get analogue expressions both for wheels and multiwheels we had to use for wheels in place of quantity q quantity s = k/2, giving fractional numbers for odd wheels. It would be interesting to ask then what pre-wheel stands behind w 1/2 . It might be multiedge or edge adjacent to loop.

Both classes, wheels and multiwheels have m = 2n -2 edges. Let us notice that to the class of graphs G n,2n-2 belong these quadrangulating projective plane. In next section we show how this fact turns crucial for multiwheels in generalizing them for projective plane.

Another question would be how to designate indices in multiwheel w q1,...,q k if we wanted to take into account type of sections standing behind corresponding indices. We have four types of sections, therefore we have to equip this index with this additional information. One way would be to use four colors for indices. Other way would be to supply index with diacritic sign, say, w 1-3-3-3-3 or w 1 3333 for w 13333 with types of sections in augmented order.

Grötzsch graph, Mycielski's Construction and 4-critical projective plane multiwheels

We start with observation that Grötzsch graph [START_REF] Bondy | Graph Theory[END_REF] may be considered as edge sum modulo two of five wheels w 1 and one wheel w 2 , see fig. Grötzsch graph is 4-critical and it quadrangulates projective plane, see fig. 7. Indeed, it has 11 vertices and 20 edges, i.e., it belongs to 2n -2-edges-class of graphs, and fig. 7 shows how this embedding on projective plane is performed. If in place of Grötzsch graph formed as 5w 1 + w 2 we take only three plus one wheel, we get graph that is isomorphic to the base graph w 13 . Taking this fact into account, we designate this graph g 111 or g 13 and traditional Grötzsch graph as g 11111 or g 15 . First, let us notice that both graphs are 4-critical, both quadrangulate projective plane, first being planar, bet second -projective planar. We might ask -are all graphs qw 1 + w q belonging summed according multiwheel summation pattern 4-critical? The answer is quite obviously positive, and we express the fact in the lemma what follows. We say that sum qw 1 +w q , k = 2q+1, q > 0, modulo two is got according multiwheel pattern if k wheels w 1 each looses two edges and w q looses k edges. This class of graphs we call Grötzsch class.

Lemma 6 For k = 2q + 1, q > 0, resulting graphs from kw 1 + w q summed according multiwheel pattern are 4-critical and quadrangulate projective plane.

Of course, for k = 1, 2, we get the base graph and Grötzsch graph, which are 4-critical, and further graphs are 4-critical due to symmetry.

This class kw 1 + w q , extending the base graph and the Grötzsch graph may be got as the Mycielski's Construction [START_REF] Bondy | Graph theory with Applications[END_REF], page 130. For that we are to take 3critical graph, i.e., arbitrary odd cycle C k , and apply Mycielski's Construction. Thus we see that if in Mycielski's Construction we replace each new got k-critical graph with arbitrary k-critical graph then we should receive k + 1-critical graph, which fact follows from the proof of the Mycielski's Construction's applicability to get k-critical graphs, see [START_REF] Bondy | Graph Theory[END_REF]. Natural question would arise does there exist Mycielski's Construction's generalization that works backwards too, i.e., that each k + 1-critical graph has as antecedent k-critical graphs in terms of this or similar construction. In order to include in Mycielski's Construction previous planar class we are to allow to match previous graphs m edges with m new vertices plus extra vertex. This would work for step from 3-critical to 4-critical graphs, and give just our plane class of multiwheels.

Further we are going to build more multiwheels, but the previous class should be the only that were quadrangulating projective plane.

Further we generalize projective planar multiwheels similarly as in case plane multiwheels, i.e., sections of w 1 may be replaced with arbitrary odd wheels. Fig. 8 shows simplest properly projective plane multiwheel q 112 . Construction 7 Let us take odd in number (k = 2q + 1) odd wheels and one wheel w q . Let us take in each of first wheels two proximal spikes and rim edge so that they do not form triangle, and middle spike edge match with central wheel w q , and other two chosen edges (spike and rim edge) match in cyclical sequence of wheels.

The resulting graph built according construction 7 belongs to 2n -2 edges class and is 4-critical. We call the resulting graph multiwheel similarly to those planar ones. Proof. Let us use the fact that the base graph belongs to Grötzsch class, and the construction for the base graph extended with non-planar section (see fig. 5) may be used for Grötzsch class in the whole.

Let us end this section with one more theorem.

Theorem 9 Multiwheel quadrangulates projective plane only if it belongs to Grötzsch class.

Proof. The only subclass to be considered is plane multiwheels with simple sections, excluding the base graph, i.e., w 1q , q > 1. It suffices to consider the minimal graph from the class w 12 . For graph to quadrangulate surface it is necessary that every edge goes into at least two square cycles. But the edge of w 15 that is incident to vertices of degree three and four doesn't fulfil this condition.

Octahedral theorem

Let us formulate what we call octahedral theorem for the plane multiwheels.

As was told in introduction, minor bracket works for the base graph, i.e., < O -, O; w 13 > is true: O -≺ w 13 , O ≺ w 13 and O -≺ O. It easily follows from facts that the base graph is only vertex 3-connected, i.e., it has triples of separating vertices, as long as octahedron graph doesn't have. This argument directly applies to plane multiwheels in general, because they are built allowing triples of separating vertices for each section, that excludes possibility for O to be minor. Both plane and projective plane multiwheels have the base graph as their minor. Besides, Grötzsch graph doesn't have O as minor. Indeed, it has 5 cubic vertices, which may be separated with triple of vertices, and adjacent to central hub vertex, and remaining 5 vertices aren't sufficient to hold O as minor. This argument easily generalizes to Grötzsch class in the whole. It only remains to persuade oneself that it works for projective plane multiwheel in general. And again, sections that are differing from simple ones can be separated by triples of vertices, see fig. 8,b. Thus, we have done with the proof.

Let us formulate the fact for arbitrary multiwheels as theorem.

Theorem 10 Minor bracket < O -, O > works for both plane and projective plane multiwheel graph classes.

Fig. 1 .

 1 Fig. 1. In the graph w1 3 (left) contracting thick edge we get octahedron graph without an edge (O -, right). Adding dotted edge we get octahedron graph O.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Forming graph w1 3 (b) by summation modulo two of edges of three graphs W3 (a): 3W3 = w1 3 . Edges that annihilate (one in one pair) by summation modulo two are depicted as dotted (b). In c) w1 3 is received by operation O + W3 where triangle (depicted dotted) is annihilated in both graphs.

.Fig. 6 .

 6 Fig.6. Example of non-planar graph that is 4-critical: Grötzsch graph. It may be built as edge sum 5W3 + W5 modulo two. Annihilated edges are depicted as dotted for one section standing for W3 and for W5.

Fig. 7 .

 7 Fig. 7. Grötzsch graph on projective plane. It quadrangulates projective plane. Compare [6].

Fig. 8 .Theorem 8

 88 Fig. 8. a) Graphs w111 and g111 are isomorphic; b)simplest properly projective plane multiwheel with minimal edge number q112.

  lower. Clearly, w 13 belongs to class of graphs with 2n -2 edges by n vertices, where all graphs quadrangulating projective plane should belong. Besides, w 13 is selfdual, if considered plane, but of course not

	a)	b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	. . . . . . . . . . . . . . . . . . . . . . . . . . .......	c)

  1-1-1 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	. . . . . . . . . . . . . . . . . . . . . . . . . . .	. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	Fig.