B Gupta

Y Hebrard

J B Coudert

B Ducharne

Characterization and local phenomenological model of Barkhausen noise under mechanical and magnetic excitations INTERMAG 2018 -GT 03 -27/04/2018 -Hysteresis Modelling I

Keywords: Barkhausen noise, magnetic hysteresis loop, model, hysteresis, mechanical stress

Post-treatment and rescaling, it is possible to plot local hysteresis cycles from the measurement of local magnetic Barkhausen noise. If the material is homogeneous and if similar excitation conditions are imposed, the local hysteresis cycles obtained are comparable to the classical magnetic hysteresis cycles B(H) (the cross-section magnetic average induction B as a function of the surface tangential excitation field H). These local Barkhausen noise hysteresis cycles give interesting clues about the evolution of the microstructure of the magnetic material (internal stresses, level of degradation etc.). This makes it an indispensable tool for the nondestructive evaluation of ferromagnetic steels. In this paper, a phenomenological modeling of the Barkhausen noise from the local modeling of B subjected to an excitation field H and/or a uni-axial mechanical stress T is proposed. The objective is to provide an absolute quantification of the internal residual stresses because of the Barkhausen noise measurements.

MAIN TEXT

The use of non-destructive micro-magnetic control techniques such as the measurement of magnetic Barkhausen noise has recently increased exponentially in the industrial environment. This renewed interest is mainly due to the improvement of analog and digital signal processing techniques that have allowed this quality control type to be integrated into the production lines in real time [START_REF] Theiner | 3MA-Testing Equipment, Application, Possibilities and Experiences. Characterization of materials[END_REF]. Micromagnetic techniques can be used to trace interesting properties of the tested samples (hardness, residual stresses, grinding burns etc.) [START_REF] Gupta | A space discretized ferromagnetic model for non-destructive eddy current evaluation[END_REF]. Barkhausen noise owes its origin to the magnetization processes of a ferromagnetic material. It is established that a ferromagnetic material, even in a demagnetized state is magnetically divided into finite regions called as domains. Each domain is characterized by its own direction and orientation of magnetization. The process of magnetization consists of converting this multi-domain state into a single domain characterized by a direction and a sense of magnetization very close to the external magnetic field H imposed on the material. This process is not continuous but consists of small discrete variations: jump from a pining defect (micro-structural obstacles like precipitation) to another defect of the domain walls (domains boundaries) [START_REF] Krause | Anisotropic flux density dependence of Magnetic Barkhausen Noise in oriented 3% Si-Fe steel Laminates[END_REF]. This local variations of magnetization states are known as "Barkhausen event". These changes induce pulsed eddy currents near the moving domain walls that develop in all space directions. These abrupt changes in the magnetic structure through the material also induces local and rapid flux variations that can be easily measured using a dedicated local micro magnetic field sensor. In recent publications dealing with the topic, many authors define a new quantity called Magnetic Barkhausen energy (MBNenergy) [START_REF] Ducharne | Characterization and modeling of magnetic domain wall dynamics using reconstituted hysteresis loops from Barkhausen noise[END_REF]. The MBNenergy is obtained from the temporal integration of the square of the local magnetic field sensor voltage output (equation 1). This "post-treatment" technique can be used to trace a hysteresis cycle MBNenergy(H) from the local measurement of magnetic Barkhausen noise. At low frequency, if the material is homogeneous and isotropic and after rescaling thanks to a coefficient ν, the Barkhausen cycles obtained are very similar to the usual hysteresis cycles B(H).

( )

2 0 ( ) . . ( 1
) T energy Barkhausen dH MBN H sign V dt dt ν   =     
These observations confirm the microscopic origin (wall movements) of the magnetic hysteresis under low frequency excitation. In the extended version of the article, the numerical scheme will be described that is used to return the temporal Barkhausen noise envelope under both magnetic and uniaxial mechanical stress. This numerical scheme is based on both local magnetic hysteresis model and an inverse method linked to magnetic Barkhausen energy.

Figure 1 .

 1 Figure 1. Magnetic Barkhausen noise (left) and simulated Barkhausen envelop under both magnetic and uniaxial mechanical stress excitations.

The good simulation results displayed in Fig. 1 are very promising. By correctly anticipating in simulation the material behavior, we will have access to some microstructural information (such as precipitations, local residual stresses). From a non-destructive evaluation point of view, this set of information will allow us to anticipate some failures and ageing issues.


Keywords: adjustment for simplicity, empirical Bayes methods, p value calibration, reproducibility crisis, strength of statistical evidence

Occam's razor suggests assigning more prior probability to a hypothesis corresponding to a simpler distribution of data than to a hypothesis with a more complex distribution of data, other things equal. An idealization of Occam's razor in terms of the entropy of the data distributions tends to favor the null hypothesis over the alternative hypothesis. As a result, lower p values are needed to attain the same level of evidence.

A recently debated argument for lowering the significance level to 0.005 as the p value threshold for a new discovery and to 0.05 for a suggestive result would then support further lowering them to 0.001 and 0.01, respectively.

Introduction

In their arguments to redefine statistical significance from 0.05 to 0.005, [START_REF] Benjamin | Redefine statistical significance[END_REF] added to the "large and ever-increasing literature on the use and misuse of significance tests," with much of it centering on the fact that the p value is not defined as a probability that the null hypothesis is true (Cox, 2006, pp. 41-42), and the resulting difficulty in interpretation [START_REF] Schervish | P values: What they are and what they are not[END_REF][START_REF] Royall | Statistical Evidence: A Likelihood Paradigm[END_REF][START_REF] Efron | Scales of evidence for model selection: Fisher versus Jeffreys[END_REF][START_REF] Goodman | Commentary: The p-value, devalued[END_REF]. The fact that p values can be automatically calculated makes them ubiquitous in scientific reports but also contributes to the perception that they are more objective than posterior probabilities.

Rather than choosing between significance testing and Bayesian methods, many researchers provided ways to calibrate a p value for interpretation as a Bayes factor, enabling interpretation as a posterior probability of null hypothesis truth when a prior probability can be specified. Such calibrations support not only the intuition of [START_REF] Fraser | Inference for bounded parameters[END_REF] and others that an extremely low p value indicates strong evidence against the null hypothesis but also support arguments against considering p values near 0.05 as indicative of strong evidence (e.g., [START_REF] Sellke | Calibration of p values for testing precise null hypotheses[END_REF]. See [START_REF] Held | On p-values and bayes factors[END_REF] for a recent review.

Under H 0 , the null hypothesis, a p value p (X) is a random variable with the uniform probability density function f 0 , where X is the random vector that models the observable sample. With x as the observed sample, p (x) then denotes the observed p value. Under the alternative hypothesis H 1 , the probability density function of p (X) is denoted by f 1 . Thus, p (X) ∼ f = π 0 f 0 +π 1 f 1 , where π 0 = P (H 0 ) and π 1 = P (H 1 ) are the prior probabilities of the null and alternative hypotheses, respectively (π 0 + π 1 = 1). B (p (x)) = f 0 (p (x)) /f 1 (p (x)) is called the Bayes factor in favor of the null hypothesis over the alternative hypothesis since 1 multiplying it by the prior odds yields the posterior odds, that is,

P (H 0 |p (x)) P (H 1 |p (x)) = B (p (x)) π 0 π 1 , (1)
by Bayes's theorem,

P (H 0 |p (x)) = π 0 f 0 (p (x)) f (p (x)) . (2)
Guided by Occam's razor, Section 2 proposes a method of adjusting the prior probability π 0 for the simplicity of H 0 relative to that of H 1 . That is accomplished indirectly by replacing an estimate of the Bayes factor such as

B ≈ 1.65 |z (x)| exp -z 2 (x)
2 with a simplicityadjusted estimate such as

|z (x)| B (p (x)) ≈ 1.65z 2 (x) e -z 2 (x) 2 , (3)
where z (x) is a z-score corresponding to p (x). For example, z (x) could be the standard normal quantile of a one-sided p value testing whether a real-valued quantity is 0, or z 2 (x) could be the χ 2 1 -quantile of a p value testing whether a non-negative quantity is 0, χ 2 1 being the χ 2 distribution with 1 degree of freedom. Equation (3) may hold when |z (x)| > 1, indicating that simplicity considerations increase the Bayes factor and thus the evidence for the null hypothesis in proportion to |z (x)|, the extent that the p value is small. As a result, the p value has to be smaller when considering simplicity than otherwise to achieve a given level of evidence as quantified by a threshold of the Bayes factor. That is the basis of Section 3's case for redefining the p value threshold for statistical significance to 0.001 instead of the 0.005 level that [START_REF] Benjamin | Redefine statistical significance[END_REF] advocate.

Rather than adjusting π 0 for the simplicity of H 0 relative to that of H 1 , Section 4 adjusts a hyperprior distribution for the simplicity of the prior distribution, leading to similar results.

The rationale for simplicity adjustments is clarified in the discussion, Section 5. The methods adjusted for hypothesis simplicity are extended in Appendix A to confidence-based methods of propagating uncertainty and in Appendix B to the case of simultaneously testing multiple null hypotheses.

2 Bayes factors adjusted for the simplicity of each hypothesis 2.1 How to adjust a Bayes factor for hypothesis simplicity

Likelihood functions adjusted for hypothesis simplicity

This subsection presents the general method of adjusting a likelihood function of a parameter for the simplicity of the distributions corresponding to the parameter values. [START_REF] Bickel | Computable priors sharpened into Occam's razors, working paper[END_REF] proposed it as an aid to specifying prior distributions, quantifying simplicity in terms of Kolmogorov complexity.

If g θ (y) is the probability density of an observed sample y given a parameter value θ, then g θ (y) as a function of θ for a given y is called the unsharpened likelihood function. It is a likelihood function in the usual sense, and "unsharpened" indicates that it has not been sharpened, that is, adjusted to account for simplicity as per Occam's razor. Let g

(1) θ (y 1 ) denote the probability density of y 1 , a single value observed as a component of y. For each value of θ, the probability density function g θ is a hypothesis about the distribution of Y , a random sample.

g θ (y) (4)
as a function of θ.

The sharpness κ controls the degree to which the unsharpened likelihood function is adjusted for hypothesis simplicity. At one extreme, L 0 is the unsharpened likelihood function, ignoring hypothesis simplicity, whereas at the other, L ∞ emphasizes simplicity to the extent of ignoring y. L 1 is the natural default, and L 2 is intermediate to L 1 and L ∞ inasmuch as 2 is the harmonic mean of 1 and ∞.

Bickel (2016) and [START_REF] Bickel | Confidence intervals, significance values, maximum likelihood estimates, etc. sharpened into Occam's razors[END_REF] originally considered κ = 1 and κ ≥ 0, respectively. The case of κ < 0 is allowed as a way to implement Alexandre Patriota's suggestion to give more rather than less weight to distributions of higher entropy (private communication, May 24, 2018). Except when otherwise specified, κ = 1.

Bayes factors adjusted for hypothesis simplicity

The Bayes factor

B (p (x)) = f 0 (p (x)) /f 1 (p (x)), not incorporating the simplicity of f 0 or f 1 ,
is called the unsharpened Bayes factor. It may be sharpened to degree κ by setting y = p (x)

and g a = g

(1) a = f a , where a, the indicator of the truth of the alternative hypothesis, satisfies a = 0 under H 0 and a = 1 under H 1 . In general, f a (p (x)) is an integrated likelihood function of a since f 0 and f 1 are the prior predictive distributions of p (X) under H 0 and H 1 , respectively.

Thus, equation (4

) gives L κ (a) = e -κ S f (1) a f a (p (x)) as the κ-sharpened (integrated)
likelihood function of a and

B κ (p (x)) = L κ (0) L κ (1) = e -κ S (f 0 ) f 0 (p (x)) e -κ S (f 1 ) f 1 (p (x)) = e -κ∆ B (p (x)) (5)
as the κ-sharpened Bayes factor, where ∆ = S (f 0 ) -S (f 1 ). B κ (p (x)) recovers the unsharpened Bayes factor when κ = 0 or ∆ = 0.

Example:

The normal Bayes factor and its adjustment for hypothesis simplicity

For a running example of adjusting a Bayes factor for the simplicity of the null hypothesis and the alternative hypothesis, a Bayes factor based on the normal distribution will be used.

While that is presented in [START_REF] Held | How the maximal evidence of p-values against point null hypotheses depends on sample size[END_REF] as a lower bound on the Bayes factor, it is developed here instead as a maximum likelihood estimate of the Bayes factor in order to connect it with the fiducial methods of Appendix A and to clarify its relation to the empirical Bayes methods of multiple testing covered in Appendix B.

Let p one (x) denote a given one-sided p value and p (x) = 2 min (p one (x) , 1 -p one (x))

the corresponding two-sided p value. Alternatively, if p (x) is the p value for a χ 2 test or for another test of a null hypothesis at the boundary of the parameter space, then let p one (x) = 1 -p (x) /2 . In either case, z (x) denotes the p one (x)-quantile for the standard

normal distribution function Φ, that is, z (x) = Φ -1 (p one (x)).
Under H 0 , z (X) ∼ N (0, 1) since p (X) ∼ U (0, 1). When assuming z (X) ∼ N (0, σ 2 ) for some σ > 1 under H 1 , the maximum likelihood estimate of B (p (x)) is the normal Bayes factor,

B (p (x)) = φ (z (x)) φ ( σ (x) z (x)) = φ (z (x)) φ ( σ (x) z (x)) = σ (x) e -( 1-σ -2 (x) ) z 2 (x) 2 for |z (x)| ≥ 1,
where φ is the standard normal density function, and σ (x) is defined as σ's maximum likelihood estimate,

σ (x) = arg sup σ>0 φ (σz (x)) = |z (x)| .
Thus, the normal Bayes factor is

B (p (x)) = |z (x)| e -z 2 (x)-1 2 = √ e |z (x)| e -z 2 (x) 2 ≈ 1.65 |z (x)| e -z 2 (x) 2 . ( 6
)
Since B (p (x)) is not adjusted for hypothesis simplicity, it is unsharpened in the sense that it is an estimate of B (p (x)), a unsharpened Bayes factor in the terminology of Section 2.1. Equation ( 5) suggests the sharpened normal Bayes factor,

B κ (p (x)) = e -κ ∆ B (p (x)) ,
where ∆ = H (N (0, 1)) -H (N (0, σ 2 )). From Michalowicz et al. (2013, p. 127)

, ∆ = ln 1 -ln σ = -ln σ, leading to B κ (p (x)) = σ κ B (p (x)) = |z (x)| κ B (p (x)). The B 1 (p (x))
case appears in equation (3).

Comparisons to other Bayes factors

To compare the sharpened Bayes factor B 1 (p (x)) to its unsharpened counterpart B 0 (p (x)) ( §2.2) and to the following lower bounds of the Bayes factor, Figure 1 displays the quantities as functions of p (x). There, the universal lower bound on the Bayes factor is B (p (x)) = e -z 2 (x) [START_REF] Held | How the maximal evidence of p-values against point null hypotheses depends on sample size[END_REF].

What Figure 1 calls the conservative lower bound on the Bayes factor is Further, its proximity to B 0 (p (x)), as seen in Figure 1, also suggests that it is not too low as an estimate of the unsharpened Bayes factor.

B (p (x)) = -ep (x) ln p (x) , (7)

Posterior probability adjusted for the simplicity of each hypothesis

The posterior probability of H 0 may be calculated from its prior probability π 0 and the Bayes factor B (p (x)):

P (H 0 |p (x)) = 1 + P (H 0 |p (x)) P (H 1 |p (x)) -1 -1 = 1 + B (p (x)) π 0 1 -π 0 -1 -1 (8)
according to equations (1) and (2). Likewise, given π 0 and the κ-sharpened Bayes factor B κ (p (x)), the κ-sharpened posterior probability of H 0 is

P κ (H 0 |p (x)) = 1 + e -κ∆ B (p (x)) π 0 1 -π 0 -1 -1 = 1 + B (p (x)) e -κ S (f 0 ) π 0 e -κ S (f 1 ) (1 -π 0 ) -1 -1 (9)
since equation ( 5) has B κ (p (x)) = e -κ∆ B (p (x)) and since ∆ = S (f 0 ) -S (f 1 ). The factor in the right-hand-side that is multiplied by B (p (x)) clarifies the rationale for sharpening the Bayes factor as an elegant way to apply Occam's razor to the assessment of the prior distribution, as will be seen in Section 5.

Negligible

Weak Moderate Strong Very Strong Overwhelming W 1 ≥ 0 W 1 ≥ 1 W 1 ≥ 2 W 1 ≥ 3 W 1 ≥ 5 W 1 ≥ 7 p 0.1
p 0.01 p 0.005 p 0.001 p 0.0005 p 0.00005

Table 1: Scales of evidence for H 1 over H 0 , with the first row according to intervals of the 1-sharpened weight of evidence in bits. The scale is an adaptation [START_REF] Bickel | A predictive approach to measuring the strength of statistical evidence for single and multiple comparisons[END_REF] of the classic base-10 scales of [START_REF] Jeffreys | Theory of Probability[END_REF] to W (p (x)) ≥ 3 and W (p (x)) ≥ 5, the two base-2 scales of [START_REF] Royall | Statistical Evidence: A Likelihood Paradigm[END_REF]. Here, W 1 abbreviates W 1 (p (x)), and p abbreviates p (x). Two of the p value thresholds values are in boldface since they indicate when a result is suggestive (α = 0.01) or significant (α = 0.001) according to Section 3.

Strength of statistical evidence sharpened for hypothesis simplicity

To measure the strength of statistical evidence that the alternative hypothesis is true, let

W (p (x)) = -log 2 B (p (x)), W κ (p (x)) = -log 2 B κ (p (x)), and W κ (p (x)) = -log 2 B κ (p (x)).
Adding "unsharpened" to the term from [START_REF] Good | Studies in the History of Probability and Statistics. XXXVII A. M. Turing's Statistical Work in World War II[END_REF], W (p (x)) is the unsharpened weight of evidence in the data favoring H 1 over H 0 ; analogously, W κ (p (x)) and W κ (p (x)) are the κ-sharpened weight of evidence and the κ-sharpened normal weight of evidence, respectively.

W κ (p (x)) reduces to the unsharpened weight of evidence if κ = 0 or if the hypotheses are equally simple in the distributional sense, that is, if

H (f 0 ) = H (f 1 ).
The weight of evidence is traditionally interpreted in terms of grades of evidence like those found in Table 1, with W (p (x)) in place of its W 1 . The p value thresholds shown there are based on W 1 (p (x)), the 1-sharpened normal weight of evidence.

3 Statistical significance adjusted for the simplicity of each hypothesis [START_REF] Benjamin | Redefine statistical significance[END_REF] argued that 0.005 should be the p value threshold for saying that a result is significant enough to claim a discovery and that the 0.05 threshold should instead indicate that a result is suggestive. The argument is based largely on the Bayes factors that would be attained at each of those thresholds. Those Bayes factors were not adjusted for the simplicity of each hypothesis.

Since the simplicity-adjusted Bayes factors are higher than those not adjusted, the p value thresholds needed to attain the same Bayes factor values are lower when simplicity is taken into account. Sharpening the normal Bayes factor for simplicity leads to these adjustments in the p value thresholds:

• The statistical significance threshold of 0.005 [START_REF] Benjamin | Redefine statistical significance[END_REF] corresponds to a threshold of the unsharpened Bayes factor equal to 3.5. The sharpened Bayes factor is about 3.5 when the p value is 0.001, yielding 0.001 as the simplicity-adjusted threshold for a significant result. As that is Table 1's threshold for strong evidence, that means statistical significance would only be declared when there is strong evidence against the null hypothesis.

• For a suggestive result, the p value threshold of 0.05 [START_REF] Benjamin | Redefine statistical significance[END_REF] corresponds to a threshold of the unsharpened Bayes factor equal to 1.1. The sharpened Bayes factor is about 1.1 when the p value is 0.01, yielding 0.01 as the simplicity-adjusted threshold for a suggestive result. That being Table 1's threshold for weak evidence, a result would only be called suggestive when the evidence against the null hypothesis p (x) ≤ 0.001

p (x) ≤ 0.005 p (x) ≤ 0.01 p (x) ≤ 0.05 Report Significant if κ = 1 Significant if κ = 0 Suggestive if κ = 1 Suggestive if κ = 0 κ = 0 W 0 (x) ≥ 5.4 W 0 (x) ≥ 3.5 * W 0 (x) ≥ 2.7 W 0 (x) ≥ 1.1 † κ = 1 W 1 (x) ≥ 3.7 * W 1 (x) ≥ 2.0 W 1 (x) ≥ 1.3 † W 1 (x) ≥ 0.1 κ = 2 W 2 (x) ≥ 1.9 W 2 (x) ≥ 0.5 W 2 (x) ≥ 0.0 W 2 (x) ≥ -0.9
Table 2: Weights of evidence corresponding to the 0.001 and 0.005 p value thresholds for a significant result and the 0.01 and 0.05 p value thresholds for a suggestive result. The higher of each threshold (0.005 or 0.05) is that proposed by [START_REF] Benjamin | Redefine statistical significance[END_REF], and the lower (0.001 or 0.01) is what it would take for a Bayes factor sharpened at level κ = 1 to attain approximately the same weight of evidence (3.5 ≈ 3.7 or 1.1 ≈ 1.3), as matched by * for significant results and † for suggestive results. Here, W κ (x) is the base-2 logarithm of the Bayes factor in favor of the alternative hypothesis over the null hypothesis when the degree of sharpness is κ.

is weak but not negligible.

Table 2 summarizes that method of adjusting the p value thresholds for the simplicity of the null hypothesis and alternative hypothesis using κ = 1 as the amount of the adjustment, with κ = 0 corresponding to no simplicity adjustment. (The weights of evidence for κ = 2 are also shown in the table since that is an intermediate value in the sense that 2 is the harmonic mean of 1 and ∞.)

4 Posterior probability adjusted for the simplicity of each π 0 value

In Section 2, the prior probability that the null hypothesis is true was assessed as guided by the simplicity of H 0 's prior predictive distribution compared to the simplicity of H 1 's prior predictive distribution. Using π 0 as the prior probability without simplicity considerations, the simplicity adjustment to the prior was passed to the Bayes factor for the reasons to be discussed in Section 5. In contrast with thereby treating π 0 as an epistemological probability in need of assessment with help from Occam's razor, this section treats π 0 as an unknown probability that is a property of the system studied. That means the simplicity of π 0 can guide the assessment of the its higher-level prior distribution, a hyperprior distribution.

The alternative hypothesis indicator A is defined as a random variable of Bernoulli distribution Bern (1 -π 0 ), that is, P (A = 0) = π 0 and P (A = 1) = 1 -π 0 . That A is equal to 1 if H 1 is true and that is equal to 0 if H 0 is true.

Just as Section 2.1.1 quantified the simplicity of each hypothesis in terms of its entropy of p (X), this section quantifies the simplicity of each value of π 0 in terms of its entropy of A. The general method of Section 2.1.1 applies to that entropy and hierarchical model with Y = A, θ = π 0 , and g θ = g π 0 = Bern (1 -π 0 ), as follows. Let Π denote the unsharpened hyperprior distribution of π 0 . If Π is a probability density function with respect to the Lebesgue measure or a probability mass function, then the κ-sharpened hyperprior probability density or the κ-sharpened hyperprior probability mass of π 0 is Π κ (π 0 ) ∝ e -κ S (gπ 0 ) Π (π 0 ), that is,

Π κ (π 0 ) = e -κ S (gπ 0 ) Π (π 0 ) e -κ S g π 0 Π (π 0 ) dπ 0 or Π κ (π 0 ) = e -κ S (gπ 0 ) Π (π 0 ) π 0 e -κ S g π 0 Π (π 0 ) (10)
respectively. The entropy of Bern (1 -π 0 ) is

S (g π 0 ) = (-π 0 ln π 0 -(1 -π 0 ) ln (1 -π 0 )) = -ln π π 0 0 (1 -π 0 ) (1-π 0 ) ; ∴ Π κ (π 0 ) ∝ e -κ S (gπ 0 ) Π (π 0 ) ∝ π π 0 0 (1 -π 0 ) (1-π 0 ) .
The corresponding κ-sharpened hyperposterior probability density or the κ-sharpened hy-perposterior probability mass is then given by

Π κ (π 0 |p (x)) ∝ Π κ (π 0 ) (π 0 f 0 (p (x)) + (1 -π 0 ) f 1 (p (x))) .
The κ-sharpened hyperposterior probability probability that H 0 is true is, as per equation ( 8),

P Πκ (H 0 |p (x)) = π 0 Π κ (π 0 |p (x)) π 0 f 0 (p (x)) π 0 f 0 (p (x)) + (1 -π 0 ) f 1 (p (x)) (11) = π 0 Π κ (π 0 |p (x)) 1 + B (p (x)) π 0 1 -π 0 -1 -1
in the probability mass case and the analogous quantity in the probability density case.

Example. Suppose, as in Section 2.2, that f 0 = N (0, 1) and f 1 = N (0, σ 2 ), and let Π be the probability mass function on the domain {1/2, 1} of two possible values of π 0 such that Π κ (1/2) = Π κ (1) = 1/2. To compare the resulting adjustment for π 0 simplicity to the corresponding adjustment for hypothesis simplicity ( §2), Figure 2 displays two proposed sharpened posterior probability estimates and their unsharpened counterparts as functions of p (x), the two-sided p value.

Discussion

This paper draws out implications of sharpening Bayes factors and the p value thresholds and posterior probabilities that depend on such Bayes factors, where "sharpening" means adjusting them for the simplicity of distributions as motivated by Occam's razor. But why should they be adjusted for simplicity? Even if Occam's razor is applicable to distributional simplicity, why sharpen Bayes factors as opposed to prior probabilities? Figure 2: Estimated posterior probability that the null hypothesis is true as a function of the two-sided p value. The solid curves are based on π 0 = 1/2 (50% unsharpened prior probability that the null hypothesis is true); the dashed curves are based on a 50% unsharpened hyperprior probability that π 0 = 1/2 and 50% that π 0 = 1; the thicker curves represent the sharpened versions of the thinner curves. Thus, in the notation of equations ( 9) and ( 11), they are P 1 (H 0 |p (x)) (solid, thicker), P 0 (H 0 |p (x)) (solid, thinner), P Π 1 (H 0 |p (x)) (dashed, thicker), and P Π 0 (H 0 |p (x)) (dashed, thinner).

The second factor in the last term of equation ( 9) sheds light on the rationale for sharpening (adjusting for simplicity), for that factor may be interpreted as the sharpened prior odds corresponding to the κ-sharpened prior probability π 0,κ = e -κ S (f 0 ) π 0 e -κ S (f 0 ) π 0 + e -κ S (f 1 ) (1 -π 0 )

,

where κ > 0 is the degree of sharpening according to Occam's razor; -S (f 0 ) and -S (f 1 )

are the degrees of simplicity of the distribution of the p value under the null hypothesis and the alternative hypothesis, respectively. The motive for sharpening is that the simplicity of f 0 compared to that of f 1 may guide the specification of the prior probability. According to this rationale, π 0 is the unsharpened prior probability of H 0 , that is, the prior probability that would be assigned without consideration of its distributional simplicity relative to that of H 1 . Then π 0,κ is a better considered prior probability in the sense that it accounts for that relative simplicity [START_REF] Bickel | Computable priors sharpened into Occam's razors, working paper[END_REF][START_REF] Bickel | Confidence intervals, significance values, maximum likelihood estimates, etc. sharpened into Occam's razors[END_REF].

Using the sharpened prior with the unsharpened Bayes factor yields the same posterior probability as using the unsharpened prior with the sharpened Bayes factor, as may be seen in equation ( 9). Whereas the rationale says an unsharpened prior should be sharpened, arranging the mathematics for instead sharpening the Bayes factor enables much clearer comparisons with previous Bayes factors without having to consider the unsharpened prior, as in Sections 2 and 3.

That rationale applies only at a higher level to Section 4, in which π 0 is an unknown parameter of the biological system or other system studied rather than an epistemological probability. There, the epistemological probability is modeled as a hyperprior.

hypothesis to be so vague that µ = 0, making H 1 an example of a local alternative hypothesis [START_REF] Held | On p-values and bayes factors[END_REF], and σ > 1. Equation ( 6) then generalizes to

B (p (x i ) ; m) = σe -( 1-σ -2 1 ) z 2 (x) 2 ,
where σ is the maximum likelihood estimate of the standard deviation of z (x 1 ), z (x 2 ), . . . , z (x m ), with each z (x i ) related to p (x i ) in the same way as Section 2.2 related z (x) to p (x). (The case B (p (x 1 ) ; 1) then recovers B (p (x 1 )).) An empirical Bayes estimate of the local false discovery rate of the ith null hypothesis is obtained by plugging B (p (x i ) ; m) into B (p (x)) in equation ( 9).

  which[START_REF] Sellke | Calibration of p values for testing precise null hypotheses[END_REF] considered as a lower bound on B (p (x)) under p (x) ≤ 1/e and a broad condition on the hazard rate that is useful for testing simple, two-sided null hypotheses. As[START_REF] Sellke | Calibration of p values for testing precise null hypotheses[END_REF] lamented, an upper bound B (p (x)) would be more desirable than a lower bound since a sufficiently low value of the upper bound would guarantee any specified amount of evidence against H 0 , the same not being true of a lower bound. For example, B (p (x)) < 1% implies that B (p (x)) < 1%, but B (p (x)) < 1% is compatible with B (p (x)) > 1%. However, B (p (x)) is nonetheless higher than (more conservative than) B (p (x)) and the other lower bounds of Bayes factors plotted in[START_REF] Benjamin | Redefine statistical significance[END_REF].

Figure 1 :

 1 Figure 1: The estimated likelihood ratio or lower bound on the Bayes factor as a function of the two-sided p value. From highest to lowest and from darkest to lightest, the curves are the sharpened normal Bayes factor, the unsharpened normal Bayes factor, the conservative lower bound on the Bayes factor, and the universal lower bound on the Bayes factor. In symbols, they are B 1 (p (x)), B 0 (p (x)), B (p (x)), and B (p (x)), respectively.

Acknowledgments

This research was partially supported by the Natural Sciences and Engineering Research Council of Canada (RGPIN/356018-2009) and by the Faculty of Medicine of the University of Ottawa.

A Replacing maximum likelihood with confidence-based methods

A.1 A fiducial modification of the normal Bayes factor

Let σ q = |z (x)| / F -1 (q), where F -1 (q) is the qth quantile of a χ 2 1 variate. Without knowledge of the Bayes factor B (p (x)), this counterpart may be calculated for |z (x)| ≥ 1 and a q ∈ [0, 1] :

in terms of the notation of Section 2.2.

To interpret B (p (x) ; q), consider the random variable Q ∼ U (0, 1). For each q ∈ [0, 1],

the value B (p (x) ; q) is the qth quantile of B (p (x) ; Q), a random variable with a fiducial distribution that depends on p (x), the fixed p value [START_REF] Bickel | Confidence distributions applied to propagating uncertainty to inference based on estimating the local false discovery rate: A fiducial continuum from confidence sets to empirical Bayes set estimates as the number of comparisons increases[END_REF].

Thus, B (p (x)) may be estimated by B (p (x) ; Q)'s median [START_REF] Bickel | Confidence distributions applied to propagating uncertainty to inference based on estimating the local false discovery rate: A fiducial continuum from confidence sets to empirical Bayes set estimates as the number of comparisons increases[END_REF],

which equation ( 6) shows to be essentially the same as the maximum likelihood estimate.

A.2 Fiducial estimates of the posterior probability of H 0 A fiducial estimate of the posterior probability of the null hypothesis may be found by substituting B (p (x) ; 1/2) for B (p (x)) in equation ( 9). Bickel (2017) considered the unsharpened (κ = 0) versions of that and other ways to use confidence distributions and other coherent fiducial distributions to propagate the uncertainty in posterior probabilities.

Those methods are equally applicable to sharpened posterior probabilities. For example, a κ-sharpened fiducial posterior probability of H 0 is the expectation value

Appendix B bridges the gap between the empirical Bayes terminology of [START_REF] Bickel | Confidence distributions applied to propagating uncertainty to inference based on estimating the local false discovery rate: A fiducial continuum from confidence sets to empirical Bayes set estimates as the number of comparisons increases[END_REF] and the evidential terminology of this paper.

B Multiple testing and posterior probabilities as local false discovery rates

From an empirical Bayes point of view, equation ( 2)'s P (H 0 |p (x)) is an unknown local false discovery rate [START_REF] Bickel | Small-scale inference: Empirical Bayes and confidence methods for as few as a single comparison[END_REF][START_REF] Bickel | Confidence distributions applied to propagating uncertainty to inference based on estimating the local false discovery rate: A fiducial continuum from confidence sets to empirical Bayes set estimates as the number of comparisons increases[END_REF], so named because of its relation to the false discovery rate in the context of testing multiple hypotheses, replacing the p (x) with a vector of p values, one for each hypothesis tested [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF]. In empirical Bayes terminology, the P κ (H 0 |p (x)) of equation ( 9) would then be the κ-sharpened local false discovery rate.

In the setting of m simultaneous hypothesis tests, the p values are denoted by z (x 1 ), z (x 2 ), . . . , z (x m ), where each x i represents a different sample. As part of an empirical Bayes approach to hierarchical models, [START_REF] Efron | Size, power and false discovery rates[END_REF]Efron ( , 2010, pp. 72-74) , pp. 72-74) considered normal distributions of z (X) conditional on the alternative hypothesis H 1 , the proposition that the null hypothesis H 0 : z (X) ∼ N (0, 1) is false. In short, the alternative hypothesis is that H 1 : z (X) ∼ N (µ 1 , σ) with (µ 1 , σ) = (0, 1), that is, with µ 1 = 0 and/or σ = 1. In his terminology, H 0 is the "theoretical null hypothesis" since its distribution of z (X) does not depend on x, the observation.

In the example of Section 2.2, the lack of background information requires the alternative


Pablo G S Carvalho
email: pablogscarvalho@gmail.com

Philippe R B Devloo

Sônia M Gomes
email: soniag@ime.unicamp.br

On the use of balanced H(div)-L 2 pair of approximation spaces for divergence-free simulations of Stokes and coupled Stokes-Darcy equations

Keywords: Mixed Finite Elements Methods, Discontinuous Galerkin, Stokes flow, Coupled Stokes-Darcy flow, Divergence-free Simulation

The performance of different classic or more recent finite element formulations for Stokes and coupled Stokes-Darcy problems are discussed. Discontinuous, H 1 -conforming, or H(div)conforming velocity approximations can be used for Stokes flows, the formulations being presented in a unified framework. Special emphasis is given to the former H(div)-conforming formulation, for which only the tangential velocity components require a penalization treatment. For incompressible fluids, this method naturally gives exact divergence-free velocity fields, a property that few schemes can achieve. When combined with classic mixed formulation for Darcy's flows, a strongly conservative scheme is derived for the treatment of coupled Stokes-Darcy problems. Unlike other methods existing in the literature, this technique can use the same combination of approximation spaces in both flow regions, and simplifies the enforcement of the coupling Stoke-Darcy interface conditions. All the methods are implemented using an object-oriented computational environment, and typical test problems are simulated to illustrate the properties of the different approximations, verifying errors, and rates of convergence. Then, an application to the simulation of a model representing selfcompacting concrete flow around reinforcing bars is presented. An homogenization technique is applied to interpret the reinforced bar domain by a Darcy's law, while a Stokes flow is considered in the remaining domain.

Introduction

In this paper, the interest is on finite element approximations of Stokes and coupled Stokes-Darcy problems, which appear in considerable number of applications [START_REF] Discacciati | Navier-Stokes/Darcy coupling: modeling, analysis, and numerical approximation[END_REF]. For instance, Darcy's flows occur for incompressible fluids in a rigid saturated porous media, appearing in various industrial processes. The model is formed by the mass balance of the fluid, and by a constitutive equation (Darcy's law), which relates the mean velocity of the fluid within the pores with a potential field (pressure) in terms of the permeability tensor and the fluid viscosity [START_REF] Chen | Computational methods for multiphase flows in porous media[END_REF]. Stokes problems are used to model free flow situations, where the fluid velocity is very slow, the viscosity is very large, or the length-scale of the flow is very small. This model is composed of a momentum balance equation, an equation for the conservation of mass, which are augmented with constitutive and state equations [START_REF] Slattery | Advanced transport phenomena[END_REF]. They have been used to understand lubrication [START_REF] Lautrup | Physics of continuous matter: exotic and everyday phenomena in the macroscopic world[END_REF], but also appear in a wide class of other fluid flow problems.

Situations where both kinds of flows occur simultaneously, in adjacent regions, are also of interest. For instance, this is the case of the modeling of the interaction of flows in wells and reservoirs, and of flows in fractured porous media, commonly found in Petroleum Engineering, or in modelling the interaction between surface (rivers) and groundwater (aquifer). Free flow channel confined by porous walls is a feature of many of the natural and industrial settings [START_REF] Hanspal | Numerical analysis of coupled Stokes/Darcy flows in industrial filtrations[END_REF]. Self-compacting concrete flow around reinforcing bars can also be represented by a coupled Stokes-Darcy problem, where the reinforced bar domain is represented by a Darcy's law obtained by homogenization, while a Stokes flow is considered elsewhere [START_REF] Kolařík | Computational homogenization of fresh concrete flow around reinforcing bars[END_REF].

Since its early beginning, in the sixties, the history of Finite Element Methods (FEM) shows a variety of formulations for Darcy and Stokes problems, which have been proposed and studied for different types of approximation spaces constructed for the variables involved. For both problems, pressure approximations are in L 2 (Ω) (i.e. they may be discontinuous) but the natural context for Stokes velocities is of continuous functions over element interfaces (H 1 -conforming) [START_REF] Girault | Finite element methods for Navier-Stokes equations: Theory and algorithms[END_REF], and the fluxes only have continuous normal components (H(div)-conforming) for Darcy's problems [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. There is a variety of non-conforming methods adopting discontinuous velocity spaces, as in Discontinuous Galerkin or in hybrid methods for both type of flows [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF][START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF][START_REF] Cockburn | Devising HDG methods for Stokes flow: an overview[END_REF], just to cite a few. There are also several works in the literature dealing with coupled Stokes-Darcy problems, and we refer to [START_REF] Correa | A unified mixed formulation naturally coupling Stokes and Darcy flows[END_REF][START_REF] Chen | H(div)-conforming finite element methods for the coupled Stokes and Darcy problem[END_REF][START_REF] Márquez | Strong coupling of finite element methods for the Stokes-Darcy problem[END_REF] for extensive bibliography review on this topic. Recently, it has regained relevance and is attracting new attention of numerical analysts, and improved numerical methods have been proposed [START_REF] Gatica | Analysis of the HDG method for the Stokes-Darcy coupling[END_REF].

One purpose of the current work is to be an extension of the recent numerical survey in [START_REF] Forti | A comparative numerical study of different finite element formulations for 2D model elliptic problems: continuous and discontinuous Galerkin, mixed and hybrid methods[END_REF], where different schemes for Darcy's problems are considered. That comparison study includes classic H 1 -conforming Galerkin method, and mixed methods, using well known or some more recent H(div)-conforming flux approximations. Different formulations using discontinuous approximations are also studied, namely, a Discontinuous Galerkin formulation and a primal hybrid method. All the formulations have been implemented and numerically compared in terms of accuracy and computational performance. Now, in the context of FEM applied to stationary Stokes problems, several methods are discussed, using an unified representation form, where the variables are, in principle, defined in broken spaces. Then each method is stated as a particular case of the general formulation, by imposing an specific continuity constraint for the approximations of velocity fields. However, as requested in any mixed formulation, in order to have stable and accurate simulations, the adopted approximation spaces for each variable should be chosen properly (verifying the inf-sup condition). In such context, the following methods for the Stokes problem shall be implemented and compared in terms of accuracy:

• Classic H 1 -conforming methods occur when continuous approximations are considered for the velocity variable [START_REF] Girault | Finite element methods for Navier-Stokes equations: Theory and algorithms[END_REF].

• Without assuming any continuity constraint in the approximations of the velocity and pressure variables, we are in the context of discontinuous Galerkin methods. [START_REF] Toselli | hp discontinuous Galerkin approximations for the Stokes problem[END_REF][START_REF] Rivière | Analysis of a discontinuous finite element method for coupled Stokes and Darcy problems[END_REF][START_REF] Lazarov | Stabilized discontinuous finite element approximations for Stokes equations[END_REF]. For this case, extra stabilization and/or symmetry terms need to be incorporated in the general formulation.

• When only the continuity of the normal components of the velocity field is enforced over element interfaces, H(div)-conforming approximation spaces are applied, and penalization is introduced just on the velocity tangential components over element interfaces.

In this sense, this method stays in between the two previous ones concerning regularity of velocity approximations. For Stokes flows, the use of this kind of approximation has been carried out in [START_REF] Wang | New finite element methods in computational fluid dynamics by H(div) elements[END_REF][START_REF] Wang | A robust numerical method for Stokes equations based on divergence-free H(div) finite element methods[END_REF], and can be applied for any pair of compatible space configuration in H(div)-L 2 for Darcy's problems.

Special emphasis shall be given to the former H(div)-conforming formulation for Stokes flows, since the discrete velocity field will be exactly divergence-free, assuming that the fluid is incompressible. Another important feature of this approach is that it offers a more dynamic treatment of boundary conditions than the standard Galerkin methods, since the normal component of the velocity is set as an essential boundary condition and is strongly enforced, but the tangential component of the velocity is treated as a natural boundary condition and is weakly imposed. Thus, this approach empowers the method for problems with discontinuous boundary conditions. It should be observed that in the context of hybrid methods, the idea of employing H(div) elements in order to get exactly divergence-free velocity has been explored in [START_REF] Carrero | Hybridized globally divergence-free LDG methods. Part I: The Stokes problem[END_REF][START_REF] Cockburn | Divergence-conforming HDG methods for Stokes flows[END_REF][START_REF] Lehrenfeld | High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows[END_REF].

For coupled Stokes-Darcy problems, we also adopt H(div)-conforming flux approximations on both sides, a method that has been analyzed in [START_REF] Kanschat | A strongly conservative finite element method for the coupling of Stokes and Darcy flow[END_REF][START_REF] Chen | H(div)-conforming finite element methods for the coupled Stokes and Darcy problem[END_REF]. On Stokes's part, the H(div)-conforming method just described above is used, and in Darcy's part, a classic mixed formulation is applied [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. In order to fit the fluids between the two domains, experimental conditions presented in [START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF] shall be considered, relating tangential component of fluid velocity and shear stress on Stokes-Darcy interface.

As emphasized in [START_REF] Kanschat | A strongly conservative finite element method for the coupling of Stokes and Darcy flow[END_REF], the most important property of this method is its capability of solving strongly the divergence-free equation, a task that few schemes can accomplish. Furthermore, the same pair of approximation spaces for velocity and pressure representations can be used in both flow regions, a characteristic not shown by most methods designed for the simulation of coupled Stokes-Darcy problems. Another advantage is that the Beavers-Joseph-Saffman condition is easier to enforce, since the bi-linear form only involves the tangential component of the velocity over the Stokes-Darcy interface.

A object oriented programming environment, called NeoPZ [START_REF] Devloo | Pz: An object oriented environment for scientific programming[END_REF], was used for the computational implementation of all the presented methods. Developed at LabMeC (Laboratory of Computational Mechanics), at the State University of Campinas, Brazil, it can be freely downloaded from http://github.com/labmec/neopz.

The organization of the present paper starts by setting the main notation in Section 2. The formulations for the Stokes problem are defined, discussed and computationally verified in Section 3. Next, the coupled Stokes-Darcy problem is considered in Section 4. In Section 5 the methods discussed in the previous sections are used to treat a fluid governed by the Stokes equation passing through obstacles in its original form or by an reduced Stokes-Darcy homogenized model. Finally, in Section 6, we conclude the paper with some final remarks.

Notation

Let Ω ⊂ R 2 be an open polygonal domain with border ∂Ω, and unit normal n exterior to Ω. Shape-regular partitions T = {Ω e , e = 1, • • • , n el } of Ω, formed by affine quadrilateral elements, shall be considered. The set Γ formed by all element edges E is called the mesh skeleton, and Γ int = {E ∈ Γ : E ⊂ Ω} denotes the set of internal edges. To each interior edge E, once and for all, a unit normal vector n E and a tangent vector τ E are associated so that { n E , τ E } form a right-hand coordinate system. If E is a boundary edge, then n E = n is the unit normal exterior to Ω. Over interfaces E ∈ Γ int between two elements, Ω 1 and Ω 2 , jump and average operators of a function v are formally defined as

v E = v 1 | E -v 2 | E , v E = 1 2 (v 1 | E + v 2 | E ) ,
where v i = v| Ω i . For boundary edges E ⊂ ∂Ω, jump and average are function traces over E.

Local approximations shall be defined in terms of one of the following polynomial spaces:

• P k or Q k : scalar polynomials of total degree k or maximum degree k in each variable.

• Q m,n : scalar polynomials of maximum degree m in x, and n in y.

• - → P k = [P k ] 2 or - → Q k = [Q k ] 2 :
vector polynomial spaces of total degree k or maximum degree k in each component.

Vector spaces

-→ V are piece-wise defined over the elements Ω e ∈ T in terms of local polynomial approximations -→ V e . They can be of the following kinds:

- → V d = { ϕ ∈ [L 2 (Ω)] 2 ; ϕ| Ωe ∈ - → V e }, - → V c = { ϕ ∈ [H 1 (Ω)] 2 ; ϕ| Ωe ∈ - → V e }, - → V div = { ϕ ∈ H(div, Ω); ϕ| Ωe ∈ - → V e }.

It is clear that for

E ∈ Γ int ϕ E = 0, if ϕ ∈ - → V c ( ϕ • τ e ) E τ E , if ϕ ∈ - → V div , since ϕ • n E = 0 for H(div)-conforming vector fields.
The scalar approximation spaces for the pressure are also piece-wise defined in terms of local polynomial spaces Ψ e . Globally, the spaces may be continuous or discontinuous

Ψ c = {ϕ ∈ L 2 0 (Ω) ∩ H 1 (Ω); ϕ| Ωe ∈ Ψ e }, Ψ d = {ϕ ∈ L 2 0 (Ω); ϕ| Ωe ∈ Ψ e }, where L 2 0 (Ω) = {ϕ ∈ L 2 (Ω); Ω ϕ dΩ = 0}.

Stokes flows

The first goal of the the current work is to present three methods for the Stokes problem in a unified form. In principle, the variables are defined in broken spaces. Then each method is stated as a particular case of the general formulation, by imposing a specific continuity constraint for the approximations of velocity fields: full discontinuous, continuous or divergence-conforming. These methods are discussed for different approximation space configurations, and numerically compared with respect to orders of convergence for velocity and pressure variables.

Consider the problem of finding the pressure p and the fluid velocity u, such that:

-∇ • T( u, p) = f in Ω, (1)
∇ • u = 0 in Ω, (2)
u = û over ∂Ω, (3)
where

f ∈ [L 2 (Ω)] 2 , û ∈ [H 1/2 (∂Ω)] 2
, and T( u, p) = 2µ D( u) -pI, D( u) = 1 2 (∇ u + ∇ u T ) are Stokes stress and symmetric strain rate tensors, respectively, µ being the constant fluid viscosity. We assume that ∂Ω û • η ds = 0.

Weak formulations

In this section, three different weak formulations for Stokes problem shall be presented in an unified form. They differ in the type of scalar piece-wise polynomial spaces Ψ used to approximate the pressure variable, and in the type of polynomial vector spaces for the velocity fields. The Discontinuous Galerkin Method (DG) [START_REF] Toselli | hp discontinuous Galerkin approximations for the Stokes problem[END_REF][START_REF] Girault | A discontinuous Galerkin method with nonoverlapping domain decomposition for the Stokes and Navier-Stokes problems[END_REF][START_REF] Rivière | Analysis of a discontinuous finite element method for coupled Stokes and Darcy problems[END_REF] considers broken spaces of type -→ V d and Ψ d . The other two methods can be interpreted as particular DG cases, by imposing specific continuity constraints: the classic H 1 -conforming method [START_REF] Girault | Finite element methods for Navier-Stokes equations: Theory and algorithms[END_REF][START_REF] Stenberg | Mixed hp finite element methods for problems in elasticity and Stokes flow[END_REF] occurs when continuous spaces of type -→ V c are used, and H(div)-conforming formulations are obtained by applying spaces of type -→ V div , as the ones analyzed in [START_REF] Wang | New finite element methods in computational fluid dynamics by H(div) elements[END_REF]. Given any kind of these approximation spaces -→ V and Ψ, the general variational formulation for the Stokes model problem can be written in the mixed form consisting in finding a pair of functions { u, p} ∈

- → V × Ψ verifying, for all { ϕ, ϕ} ∈ - → V × Ψ, a S ( u, ϕ) + b S ( ϕ, p) = f ( ϕ) + (BC), (4) b
S ( u, ϕ) = ∂Ω ϕû • n ds, (5)
where f ( ϕ) = Ω ϕ • f dΩ, and the bi-linear forms are

b S ( u, ϕ) = - nel e=1 Ωe ∇ • u ϕ dΩ e + E∈Γ int E ϕ u • n E ds, (6)
a S ( u, ϕ) = 2µ nel e=1 Ωe D( ϕ) : D( u) dΩ e -2µ E∈Γ E D( u) n E • ϕ ds -β E D( ϕ) n E • u ds + E∈Γ γ E u • ϕ ds. (7)
The parameter β = ±1 is taken to have a non-symmetric or a symmetric bi-linear form a S , and the term penalized by γ is usually included for stability purposes. The boundary term (BC) shall be defined separately for each formulation.

Discontinuous Galerkin Method

The DG method corresponds to take in the weak variational formulation ( 4)-( 5) test and trial functions in

- → V d × Ψ d .
The boundary term takes the form

(BC) d = 2µβ ∂Ω D( ˆ ϕ) n • û ds + ∂Ω γ û • ϕ ds. (8)
In [START_REF] Toselli | hp discontinuous Galerkin approximations for the Stokes problem[END_REF], the suggestion is to use penalization of the form γ| E = γ 0 k 2 |E| , k being the degree of polynomials chosen for the velocity approximation. The parameter γ 0 > 0 has to be chosen by the user in order to get stability (inf-sup condition). It has been shown in [START_REF] Rivière | Analysis of a discontinuous finite element method for coupled Stokes and Darcy problems[END_REF][START_REF] Lazarov | Stabilized discontinuous finite element approximations for Stokes equations[END_REF] that, for the symmetric case (β = -1), γ 0 should be sufficiently large, and bounded by a positive value. For the non-symmetric case (β = 1), taking a fixed constant γ 0 > 0 is sufficient. This kind of constraint is standard for the symmetric interior penalty Galerkin (SIPG) method applied to Poisson problems, but for the non-symmetric interior penalty Galerkin (NIPG) case, γ 0 = 0 is allowed for such problems, which is not the case for Stokes problems.

H 1 -conforming formulation

For vector test and trial functions in -→ V c , the jump operator vanishes at Γ int . Furthermore, assuming u| ∂Ω = û and taking ϕ ∈ -→ V c 0 , i.e, vanishing on ∂Ω, thus the weak formulation becomes the classic H 1 -conforming scheme

a c S ( u, ϕ) + b c S ( ϕ, p) = f ( ϕ), (9) b c S ( u, ϕ) = 0, (10)
for all in ∀ { ϕ, ϕ} ∈ -→ V 0 × Ψ, for which the bi-linear forms a c S ( u, ϕ) and b c S ( u, ϕ) reduce to

a c S ( u, ϕ) = 2µ Ω D( ϕ) : D( u) dΩ, b c S ( u, ϕ) = - Ω ∇ • u ϕ dΩ.

Divergence-conforming formulation

Consider the weak formulation ( 4)-( 5) based on a vector space -→ V div and a pressure space Ψ d such that the pair -→ V div , Ψ d is compatible for mixed formulations for Poisson problems based on T . This means that their local polynomial spaces are related by ∇ • -→ V e = Ψ e . Consequently, some terms of the bi-linear forms can also be modified by using the continuity of normal components of the velocity field, expressing the jumps and average operators only in terms of tangential components, as described in [START_REF] Wang | New finite element methods in computational fluid dynamics by H(div) elements[END_REF]. Precisely, the H(div)-

conforming formulation consists in finding { u, p} ∈ - → V div × Ψ d such that the equations a div S ( u, ϕ) + b div S ( ϕ, p) = f ( ϕ) + (BC) div . ( 11
) b div S ( u, ϕ) = 0, (12)
hold for all { ϕ, ϕ} ∈ -→ V div × Ψ d , where the corresponding bi-linear forms are

a div S ( u, ϕ) = 2µ nel e=1 Ωe D( ϕ) • D( u) dΩ e -2µ E∈Γ E D( u) n E • τ E ϕ • τ E ds -β E D( ϕ) n E • τ E • u • τ E ds + E∈Γ γ E |E| E u • τ E • ϕ • τ E ds, (13)
b div S ( u, ϕ) = - Ω ∇ • u ϕ dΩ. (14)
Similarly, the boundary term becomes

(BC) div = 2µβ ∂Ω (D( ϕ) n • τ )(û • τ ) ds + ∂Ω γ(û • τ )( ϕ • τ ) ds + ∂Ω ϕû • n ds.
It should be observed that, because of the space compatibility constraint, equation ( 12) implies that the incompressible property ∇ • u = 0 holds strongly in Ω.

Error analysis for the mixed weak formulations

The error analysis considered here holds for Stokes problems applied to homogeneous media. For simplicity, the homogeneous Dirichlet condition u| ∂Ω = 0 is considered. Let u and p be the exact solutions, and u h ∈ -→ V h and p h ∈ Ψ h their approximations based on shape-regular partitions T h , with mesh size h. Norms for L 2 (Ω) and H m (Ω) are used for error estimations, denoted by || • || and || • || m , respectively. In the context of methods using discontinuous velocity approximations, mesh dependent H 1 -norms || • || 1,h are also applied.

As is well known [START_REF] Brezzi | Mixed discontinuous Galerkin methods for Darcy flow[END_REF], the convergence of mixed methods depends not only on the approximation potential of the spaces -→ V h , and Ψ h , but also on the satisfaction of a inf-sup stability condition between them. The verification of the inf-sup condition is simpler when stabilization techniques are imposed, allowing more flexibility in the choice of the approximation spaces. Without applying stabilization, e.g. when using H 1 -conforming formulation, Method Space L 2 -error rates Divergence configurations in u in p free

DG - → Q d k Q d k-1 , N-S k k No - → Q d k Q d k-1 , S k + 1 k No - → Q d RT (k) Q d k , N-S k k No - → Q d RT (k) Q d k , S k + 1 k No H 1 - → Q c k Q c k-1 k + 1 k No - → Q c k+1 Q d k-1 k + 1 k No H(div) - → Q div RT (k) Q d k , N-S k k Yes - → Q div RT (k) Q d k , S k + 1 k Yes
Table 1: Type of space configurations used for approximations of the variables u and p for two-dimensional Stokes problems defined on convex regions, and some of their main properties.

or by just doing partial stabilization, when H(div)-conforming approximations are used, the approximation spaces are required to satisfy some particular compatibility constraints.

As shall be described in the next sections, rates of convergence for both variables are dependent on each other, and they are affected by the choice of the approximation spaces. Precisely, the smallest value between the orders of approximation of the exact solutions u in -→ V h and p in Ψ h (measured in H 1 and L 2 norms, respectively) determines the rate of convergence of the approximate velocity u h and pressure p h , simultaneously, under the corresponding norms. However, the velocity error measured by the L 2 -norm can be one unit higher, provided adjoint consistent formulations are used.

Table 1 summarizes all the methods that shall be numerically compared in the next sections, and their main properties with respect to expected orders of convergence for velocity and pressure, using the L 2 -norms, and the strongly verification of divergence-free velocity.

Stabilized full discontinuous weak formulations

Without any continuity constraint, the spaces -→ V d × Ψ d are usually formed by local polynomial spaces restricted to the elements Ω e . In [START_REF] Lazarov | Stabilized discontinuous finite element approximations for Stokes equations[END_REF], a generalized error analysis for approximate velocity and pressure variables, obtained by stabilized weak formulations (4)-( 5), has been developed. No restriction on the geometry of the elements Ω e is necessary, and quite general local finite element spaces of discontinuous functions are allowed, by just requiring them to satisfy the restriction

∇Ψ e ⊂ - → V e . ( 15
)
For the symmetric case, the penalization function γ should be sufficiently large, and bounded bellow by a positive value. For the non-symmetric case, it is sufficient to take γ > 0. Under such assumptions, the inf-sup condition holds, and assuming that P l ⊂ -→ V e and P m ⊂ Ψ e , with l ≥ 1 and m ≥ 0, then

|| u -u h || 1,h + ||p -p h || ≤ C h l || u|| l+1 + h m+1 ||p|| m+1 , ( 16
)
where C is a constant, independent of h, and || • || 1,h denotes a standard mesh-dependent H 1 -norm used in DG methods. Furthermore, assuming elliptic-regularity, holding for bidimensional convex regions, and applying the symmetric weak formulation, the following improved order of convergence holds for the velocity error measured by the L 2 -norm

|| u -u h || ≤ C h l+1 || u|| l+1 + h m+2 ||p|| m+1 . (17)
These results correspond to Theorem 4.1 and Theorem 4.2 in [START_REF] Lazarov | Stabilized discontinuous finite element approximations for Stokes equations[END_REF].

For the current applications, the following two space configurations, verifying the restriction [START_REF] Gatica | Analysis of the HDG method for the Stokes-Darcy coupling[END_REF], shall be used for quadrilateral meshes, in symmetric and non-symmetric forms:

• Space configuration of type - → Q d k , Q d k-1 : discontinuous velocity space with - → V e = - →
Q k , and discontinuous pressure space with Ψ e = Q k-1 , for k ≥ 2. Note that this particular case has been previously analyzed in [START_REF] Schötzau | Mixed hp-DGFEM for incompressible flows[END_REF].

• Space configuration of type - → Q d RT (k) , Q d k , k ≥ 1:
the pressure space is discontinuous with Ψ e = Q k , and the velocity space is discontinuous with

- → V e = Q k+1,k × Q k,k+1
, for k ≥ 1, corresponding to the local spaces used in the construction of the classic Raviart-Thomas spaces for mixed formulations of elliptic problems on quadrilateral element geometry.

Methods without or with partial penalization

These are the cases of H 1 -conforming (no penalization) or H div -conforming (partial penalization of the tangential components of the velocity) formulations. Usually, their local spaces -→ V e and Ψ e , defined on Ω e , are constructed by means of standard transformations, based on the geometric diffeomorphisms F Ωe : KΩe → Ω e , which map polynomial spaces defined on the master element KΩe to -→ V e or to Ψ e . When F Ωe is affine, the functions in -→ V e and Ψ e are polynomials as well. For such methods without (or with insufficient) penalization, the inf-sup condition requires that the approximation spaces for velocity and pressure verify some special compatibility constraints.

H 1 -conforming formulation

On this context, there exist rather general techniques for verifying the inf-sup condition and for the construction of compatible space configurations

- → V c h × Ψ h satisfying it, both with continuous Ψ h = Ψ c h or discontinuous Ψ h = Ψ d h pressure spaces.
A large collection of such methods, which are known to be stable, are described in the book by Girault and Raviart [START_REF] Girault | Finite element methods for Navier-Stokes equations: Theory and algorithms[END_REF]. The analysis in [START_REF] Stenberg | Mixed hp finite element methods for problems in elasticity and Stokes flow[END_REF] also considers several stable approximation space configurations, for quadrilateral or hexahedral geometry, with discontinuous pressure spaces Ψ h = Ψ d h . For such stable methods, classic error estimations are of type

|| u -u h || 1 + ||p -p h || ≤ C inf v∈ -→ V c h || u -v|| 1 + inf ϕ∈Ψ h ||p -ϕ|| . (18)
Furthermore, assuming elliptic-regularity, enhanced L 2 -errors for the velocity holds:

|| u -u h || ≤ Ch    inf v∈ -→ V div h || u -v|| h 1 + inf ϕ∈Ψ h ||p -ϕ||    . ( 19
)
The following cases, based on affine quadrilateral meshes, shall be compared:

• Space configuration of type - → Q c k , Q c k-1 :
these are the well known Taylor-Hood spaces (k = 2), where

- → V e = - →
Q k and continuous pressure space with Ψ e = Q k-1 , which have been analyzed in [START_REF] Stenberg | Error analysis of some finite element methods for the Stokes problem[END_REF].

• Space configuration of type - → Q c k+1 , Q d k-1 , k ≥ 2:
using discontinuous pressure space with Ψ e = Q k-1 , the velocity approximation spaces should be larger than -→ Q k , as discussed in [START_REF] Stenberg | Mixed hp finite element methods for problems in elasticity and Stokes flow[END_REF]. One stable possibility is to take

- → V e = - → Q k+1 .

Divergence-conforming formulation

Assuming velocity approximation spaces

- → V div h in H(div, Ω)
, only the tangential component of the velocity needs to be penalized at element interfaces. As described before, for general DG methods, the penalization function γ to be applied to symmetric H(div)conforming formulations should also be sufficiently large, and bounded bellow by a positive value. For the non-symmetric case, it is sufficient to take γ > 0. Using this partial penalization, the inf-sup condition for the Stokes problem can be obtained by imposing an extra compatibility constraint [START_REF] Wang | New finite element methods in computational fluid dynamics by H(div) elements[END_REF]: the pair -→ V div h × Ψ d h should form compatible spaces for the mixed formulation of Poisson problems. Under these hypotheses, Theorem 5.1 and Theorem 5.2 in [START_REF] Wang | New finite element methods in computational fluid dynamics by H(div) elements[END_REF] imply error estimations of type ( 16) and ( 17), the convergence rates depending on the particular family of compatible spaces adopted. Note that similar results hold for Navier-Stokes flows [START_REF] Cockburn | A note on discontinuous Galerkin divergence-free solutions of the Navier-Stokes equations[END_REF], using discretizations in the context of local discontinuous Galerkin (LDG) methods, based on space configurations which are compatible for the Poisson problem.

The following space configuration for affine quadrilaterals meshes shall be applied:

• Space configuration of type - → Q div RT (k) , Q d k , k ≥ 1:
this is the well known Raviart-Thomas space configuration introduced in [START_REF] Raviart | A mixed finite element method for 2-nd order elliptic problems[END_REF] for mixed methods of Poisson problems. The pressure space is discontinuous with Ψ e = Q k , and the velocity space is

H(div)-conforming with - → V e = Q k+1,k × Q k,k+1 , verifying the compatibility property ∇ • - → V e = Ψ e .
For this type of spaces, projections Π h u and Λ h p over -→ V div h and Ψ h can be defined commuting the divergence operator: ∇

• - → V div h = Λ h [∇ • u]. Furthermore, || u -Π h u|| 1,h = O(h k ) and ||p -Λ h p|| = O(h k+1 ).
With the help of such projections, error estimations of type ( 16)-( 17) hold with l = m = k.

Verification tests

A simple Stokes model, with known analytic solution, is adopted in order to verify the implemented computational code for all the described weak formulations. The domain is Ω = (0, 1) × (0, 1), and with properly chosen forcing term f , the functions u = sin(2πx) cos(2πy) -cos(2πx) sin(2πy)

, p = - 1 π cos(πx) e y , solve the problem -∆ u + ∇p = f in Ω, ∇ • u = 0 in Ω, u • n = 0 on ∂Ω.
The results obtained with the methods listed in Table 1, based on square meshes Π h with spacing h = 1/N, N = 2 i , i = 2, • • • 6, shall be compared in terms of errors measured by L 2 -norm.

On the choice of the stabilization parameter γ As already mentioned, the function γ is piece-wise defined over the skeleton by γ| E = γ 0 k 2 |E| , k being the polynomial order chosen for the approximation of the velocity. To identify good choices of the parameter γ 0 concerning stability and accuracy, a numerical study has been performed for the current test problem to evaluate the behaviour of the discontinuous and the H(div)-conforming methods based on the mesh with h = 0.125. The plots in Fig. 1 show that the accuracy of u stabilizes for increasing values of γ 0 . The symmetric DG formulations using space configurations of type

- → Q d k Q d k-1
are more sensitive to the choice of this parameter, requiring much larger values of γ 0 to reach a steady accuracy, a behaviour that is more significant as k increases. The methods using space configurations of type

- → Q d RT (k) Q d k and - → Q div RT (k) Q d
k are less sensitive to the variation of γ 0 . Based on this numerical analysis, γ 0 = 20 was considered a safe choice for all methods, and shall be adopted for all simulations presented in what follows. |E| on the accuracy of u using symmetric (S) and non-symmetric (N-S) methods for space configurations

■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ ■ Q  RT (2) div Q 2 d (N-S) ◆ Q  2 d Q 1 d (N-S) ▲ Q  RT (2) d Q 2 d (N-S) ▼ Q  RT (2) div Q 2 d (S) ○ Q  2 d Q 1 d (S) □ Q  RT (2) d Q 2 d (S) 0 10 20 30 40 50
1.E-2 1.E-3 γ0 || u  -u  h|| ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ ■ Q  RT (3) div Q 3 d (N-S) ◆ Q  3 d Q 2 d (N-S) ▲ Q  RT (3) d Q 3 d (N-S) ▼ Q  RT (3) div Q 3 d (S) ○ Q  3 d Q 2 d (S) □ Q  RT (3) d Q 3 d (S)
- → Q d k Q d k-1 , - → Q d RT (k) Q d k and - → Q div RT (k) Q d k ,
on a square mesh with spacing h = 0.125, and using k = 2 (left side) and k = 3 (right side).

Confirmation of the predicted rates of convergence

The errors in u and p, obtained with all the methods listed in Table 1 are displayed in Table 2, for k = 2 and 3. The expected convergence rates of order k + 1 for u are clearly observed in all the adjoint consistent cases. For the non-symmetric methods, stronger

N -→ Q d 2 Q d 1 (N-S) -→ Q d 2 Q d 1 (S) N -→ Q d 3 Q d 2 (N-S) -→ Q d 3 Q d 2 (S) || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || 4 2.04E -2 1.25E -1 1.92E -2 6.58E -2 4 2.01E -3 3.59E -2 1.95E -3 3.22E -2 8 2.84E -3 2.58E -2 2.51E -3 1.60E -2 8
N -→ Q d RT (2) Q d 2 (N-S) -→ Q d RT (2) Q d 2 (S) N -→ Q d RT (3) Q d 3 (N-S) -→ Q d RT (3) Q d 3 (S) || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || 4 1.38E -2 5.78E -2 1.40E -2 6.99E -2 4
N -→ Q c 2 Q c 1 -→ Q c 3 Q d 1 N -→ Q c 3 Q c 2 -→ Q c 4 Q d 2 || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || 4 2.00E -2 1.78E -2 2.00E -3 1.65E -2 4
N -→ Q div RT (2) Q d 2 (N-S) -→ Q div RT (2) Q d 2 (S) N -→ Q div RT (3) Q d 3 (N-S) -→ Q div RT (3) Q d 3 (S) || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || || -→ u --→ u h || ||p -p h || 4 1.38E -2 5.88E -2 1.40E -2 6.99E -2 4 1.36E -3 6.85E -3 1.35E -3 6.80E -3 8 1.85E -3 8.20E -3 1.86E -3 9.54E -3 8
→ Q d k Q d k-1 and - → Q d RT (k) Q d k ; H 1 -conforming method of types - → Q c k Q c k-1 , and - → Q c k+1 Q d k-1
; H(div)-conforming symmetric and non-symmetric methods of type

- → Q div RT (k) Q d k .
orders occur, greater than the expected order k, specially for the spaces of type

- → Q div RT (k) Q d k , - → Q div RT (k) Q d k , and
- → Q d 3 Q d 2 ,
for which convergence rates of order k + 1 are clearly observed, one unit higher than the predicted order k for methods with lack of adjoint consistency. Probably this fact could be justified by the size of the penalization γ| E = 20k 2 |E| , similarly to the case of non-symmetric DG methods for Poisson problems [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF], where the theoretically predicted loss of optimality in L 2 for the velocity can be compensated by imposing strong penalty. It is also observed that the choice of the approximation space for the velocity (continuous, discontinuous or divergence-conforming) does not have a significant impact on the velocity accuracy. Excepting the non-symmetric

- → Q d 2 Q d 1 case
, all the other velocity space configurations give comparable velocity accuracy for a given polynomial degree. The difference in error is smaller when using higher order schemes (e.g. k = 3). Concerning the accuracy of the pressure, the predicted order k clearly holds for both stable configurations of type

- → Q c k Q c k-1 and - → Q c k+1 Q d k-1
when using H 1 -conforming approximation spaces, and for the DG methods with spaces of type

- → Q d k Q d k-1 .
However, super-convergence of order k + 1 is observed for all methods using configurations of type

- → Q d RT (k) Q d k and - → Q div RT (k) Q d k ,
a fact that the theoretical analysis cannot explain. Similar faster convergence results for pressure have also been numerically observed in [START_REF] Cockburn | A note on discontinuous Galerkin divergence-free solutions of the Navier-Stokes equations[END_REF] for a DG method applied to a Navier-Stokes problem using

- → Q div RT (k) Q d k space configuration.
With respect to symmetric versus non-symmetric schemes, their results seem to be quite similar, except for the

- → Q d 2 Q d 1 case
, for which the symmetric case gives better pressure accuracy.

These results are represented graphically by the plots in Fig. 2. The L 2 -norms of ∇ • u are displayed in Table 3 illustrating the strongly enforced divergence-free property obtained with the

- → Q div RT (k) Q d k space configurations.
- → Q d k Q d k-1 ; H 1 -conforming method of type - → Q c k Q c k-1
; H(div)-conforming symmetric and non-symmetric methods of type

- → Q d RT (k) Q d k .
Error curves obtained with the H 1 -conforming method are shown in Fig. 3 for the velocity u (left side) and for p (right side), with stable space configurations of type

- → Q c k Q c k-1 , and - → Q c k+1 Q d k-1
, using k = 2 (top side) and k = 3 (bottom side). Since the rates of convergence of both variables are dependent on each other, the enlargement of the velocity approximation space in the -→ Q c k Q c k-1 space configuration does not imply in a better order of convergence, but the velocity errors are more than one (or even two) orders of magnitude less than the ones obtained with the Taylor-Hood method. The results for its counterpart space configuration

- → Q c k Q d k-1
, using discontinuous pressure approximations, are also include in Fig. 3. Note that the error analysis in [START_REF] Stenberg | Mixed hp finite element methods for problems in elasticity and Stokes flow[END_REF] does not apply to this framework. In fact, according to it (see Method 6 in [START_REF] Stenberg | Mixed hp finite element methods for problems in elasticity and Stokes flow[END_REF]), taking the velocity approximations in -→ Q c k , a maximal stable space combination requires discontinuous pressure with local polynomials in Q k-2 ∪ P k-1 . Despite this case, for this particular test problem, the results for

- → Q c k Q d k-1
appear to be quite stable, the velocity being in close agreement with its counterpart

- → Q c k Q c k-1 .
For k = 3, the pressure approximations of these two methods are also similar.

• • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ • Q 2 c Q 1 c ■ Q 2 c Q 1 d ◆ Q 3 c Q 1 d 2.E-2 5.E-2 1.E-1 2.E-1 1.E-2
1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

h ||u  -u  h || 3 1 3 1 • • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ • Q 2 c Q 1 c ■ Q 2 c Q 1 d ◆ Q 3 c Q 1 d 2.E-2 5.E-2 1.E-1 2.E-1 1.E-2 1.E-3 1.E-4 h ||p  -p  h || 2 1 • • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ • Q 3 c Q 2 c ■ Q 3 c Q 2 d ◆ Q 4 c Q 2 d 2.E-2 5.E-2 1.E-1 2.E-1 1.E-3 1.E-4
1.E-5

1.E-6

1.E-7

1.E-8

1.E-9

1.E-10

h ||u  -u  h || 4 1 4 1 • • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ • Q 3 c Q 2 c ■ Q 3 c Q 2 d ◆ Q 4 c Q 2 d 2.E-2 5.E-2 1.E-1 2.E-1 1.E-2
1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

h ||p  -p  h || 3
→ Q c k Q c k-1 , - → Q c k+1 Q d k-1
, and space configuration

- → Q c k Q d k-1
, based on uniform square meshes with spacing h, for k = 2 (top side) and k = 3 (bottom side).

N -→ Q d 2 Q d 1 (N-S) -→ Q d 2 Q d 1 (S) N -→ Q d RT (2) Q d 2 (N-S) -→ Q d RT (2) Q d 2 (S) ||∇ • ( -→ u --→ u h )|| ||∇ • ( -→ u --→ u h )|| ||∇ • ( -→ u --→ u h )|| ||∇ • ( -→ u --→ u h )|| 4 3.66E -1 3.69E -
N -→ Q c 2 Q c 1 -→ Q c 3 Q d 1 N -→ Q div RT (2) Q d 2 (N-S) -→ Q div RT (2) Q d 2 (S) ||∇ • ( -→ u --→ u h )|| ||∇ • ( -→ u --→ u h )|| ||∇ • ( -→ u --→ u h )|| ||∇ • ( -→ u --→ u h )|| 4
- → Q c k Q c k-1 and - → Q c k+1 Q d k-1
, the symmetric (S) and non-symmetric (N-S) DG method of types

- → Q d k Q d k-1
, and the H div -conforming method of type

- → Q div RT (k) Q d 2 , for k = 1, 2
, and square meshes with spacing h = 1/N .

Coupled Stokes-Darcy model

We present in this section a finite element formulation for the coupled model considering a Darcy's flow in a region Ω D and a Stokes flow in a region Ω S . Both regions form the computational domain Ω = Ω D ∪ Ω S , and are assumed to be polygonal, sharing the interface Γ SD = ∂Ω D ∩ ∂Ω S . Let n SD be unit normal vectors to the edges in Γ SD , pointing from Ω S to Ω D , and let τ SD be the associated positive oriented tangent vectors. For convenience, restrictions of functions to each of the flow domains are denoted by v S = v| Ω S , and v D = v| Ω S .

The coupled Stokes-Darcy model problem consists in finding u and p such that:

-∇ • T( u, p) = f in Ω S , µK -1 u + ∇p = 0 em Ω D , ∇. u = f in Ω,
where K is the permeability tensor, f ∈ [L 2 (Ω S )] 2 , and f ∈ L 2 (Ω), with f | Ω S = 0. The boundary conditions for this coupled problem are:

u = 0 in ∂Ω S \ Γ SD , u . n = 0 in ∂Ω D \ Γ SD .
Furthermore, the following conditions should be enforced at the interface Γ SD :

• Flux continuity: u S • n SD = u D • n.
• Balance of normal forces:

p S -2µ [D( u S ) n SD ] • n SD = p D .
• Beavers-Joseph-Saffman (BJS) condition [START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF]:

u S • τ SD = -2G [D( u S ) n SD ] • n SD ,
where G > 0 is an empirical coefficient.

[13] for vector spaces of type -→ V div , where general H(div)-L 2 balanced spaces are allowed, in different geometry, and numerical simulation results are presented for quadrilateral meshes, demonstrating the main qualities of the method. A recent study [START_REF] Márquez | Strong coupling of finite element methods for the Stokes-Darcy problem[END_REF] considers different forms for the combined mixed method using space configuration of type -→ V c,div for the velocity fields, combined with appropriate pressure spaces Ψ, which should be of type Ψ d on the Darcy's part, but can be of type Ψ c or Ψ d on the Stokes part, depending on the kind of H 1 -conforming space used there.

As already mentioned, the special interest on the method using balanced H(div) -L 2 space configurations for the simulation of Stokes and coupled Stokes-Darcy mixed methods comes mainly from its important capability of solving strongly the divergence-free equation, a task that few schemes can accomplish. Furthermore, for coupled flows, the same pair of approximation spaces can be used in both Stokes and Darcy's regions, a convenient characteristic not shown by most methods designed for such problems. Another advantage is that the Beavers-Joseph-Saffman condition is easier to enforce, since the bi-linear form only involves the tangential component of the velocity over the Stokes-Darcy interface.

The next example verifies the implemented methodology of these divergence-conforming methods by applying them to a Stokes-Darcy problem with known exact solution.

Verification test

The following problem was proposed by [START_REF] Chen | Weak galerkin method for the coupled Darcy-Stokes flow[END_REF]. The region is a rectangle where the Stokes domain is Ω S = (0.0, π) × ([0.0, 1.0), and the Darcy domain is Ω D = (0.0, π) × (-1.0, 0.0), with interface Γ SD = {0 < x < π, y = 0}. Viscosity coefficient µ = 1, constant permeability K = I, and coefficient for the BJS interface condition µ * = 1 are adopted. For this problem, the exact velocity and pressure fields are: The forcing functions g and f are obtained from these exact solutions. Approximations u h ∈ -→ V h , p h ∈ Ψ h are computed using uniform rectangular meshes T h with mesh sizes h xi = π/N , and

u S = v (
h y i = 2/N , for N = 2 j , j = 2, • • • , 6, and H(div) -L 2 balanced spaces - → V div Ψ d of type Q div RT (k) Q d k , k = 1, 2 and
3. For this verification test, symmetric formulation with stabilization parameter γ 0 = 12 is considered on the Stokes part. The results are represented graphically in Fig. 4, where the L 2 -rates of order k + 1 for velocity and pressure, predicted by the Raviart-Thomas mixed scheme on the Darcy's domain, have been achieved. The results for the Stokes flow are consistent with the ones obtained in the previous section, and correspoding to the predicted order k + 1 for velocity, and rates close to order k + 1 for pressure, which is higher than the predicted order k.

• • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ • Q RT (1) div Q 1 d ■ Q RT (2) div Q 2 d ◆ Q RT (3) div Q 3 d 2.E-2 5.E-2 1.E-1 2.E-1 1.E-1 1.E-2
1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

h ||u  S -u  S h || 4 1 3 1 2 1 • • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ • Q RT (1) div Q 1 d ■ Q RT (2) div Q 2 d ◆ Q RT (3) div Q 3 d 2.E-2 5.E-2 1.E-1 2.E-1 1.E+0 1.E-1 1.E-2 1.E-3 1.E-4 1.E-5 h ||p  S -p  S h || 4 1 3 1 2 1 • • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ • Q RT (1) div Q 1 d ■ Q RT (2) div Q 2 d ◆ Q RT (3) div Q 3 d 2.E-2 5.E-2 1.E-1 2.E-1 1.E-1 1.E-2
1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

1.E-8

1.E-9

h ||u  D -u  D h || 4 1 3 1 2 1 • • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ • Q RT (1) div Q 1 d ■ Q RT (2) div Q 2 d ◆ Q RT (3) div Q 3 d 2.E-2 5.E-2 1.E-1 2.E-1 1.E-1 1.E-2
1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

1.E-8

1.E-9

h ||p  D -p  D h || 4

Application

This section is dedicated to the application of the methods discussed in the previous sections to treat a Stokes fluid passing through obstacles. This kind of problem is of interest in Civil Engineering applications to represent self-compacting concrete flow around reinforcing bars [START_REF] Sandström | A two-scale finite element formulation of Stokes flow in porous media[END_REF][START_REF] Kolařík | Computational homogenization of fresh concrete flow around reinforcing bars[END_REF]. The idea is to replace the original multiscale problem, that would require strongly refined meshes to represent the small scales, by a Darcy's type flow with an homogenized permeability property, which is simpler to resolve. The procedure involves an expensive first step of Stokes simulation in order to estimate the homogeneous permeability coefficient. Then different scenarios can be solved by a reduced coupled Stokes-Darcy model, using less refined space configurations. This strategy is illustrated in this section using the following test example. In order to obtain a permeability coefficient for Darcy domain, compatible with the obstacle region of the first problem, the pressure difference ∆p num = 188, 44 between left and right extremity points of Ω 2 in the simulation of Step 1 is computed. It is shown that pressure has approximately linear variation through this second domain. Then the permeability coefficient κ num = 0, 0212 can be obtained from numerical results. For verification, this is compared with an estimated coefficient κ est by the homogenization procedure suggested in [START_REF] Tamayol | Analytical determination of viscous permeability of fibrous porous media[END_REF], where different curves are plotted for the parameters: ε = 1 -πd 2 4S 2 and κ * = κest d 2 , related to geometrical properties of the obstacles (S, surface and d, diameter). In [START_REF] Tamayol | Analytical determination of viscous permeability of fibrous porous media[END_REF], for obstacles of radius r.

Based on the comparative study shown in Table 4, the permeability value κ num can be considered and Step 2 is performed. For comparison, the pressure fields obtained in Step 1 and Step 2 are illustrated in Fig. 7. Pressure profiles at y = 2 are plotted in Fig. 8. Observe that the pressure is constant in Ω 1 (free flow domain) and then, when the flow reaches the obstacle domain, the pressure varies close to linear profile verified by the reduced coupled Stokes-Darcy model of Step 2.

Conclusions

Approximating the Stokes equations using discontinuous Galerkin methods yields consistent results with optimal convergent rates. The use of discontinuities between elements simplifies the stabilization of numerical approximations, by means of penalization, giving flexibility on the choice of approximation spaces and element geometry. They are also appealing for problems where there is imposition of discontinuous boundary conditions, and in the treatment of interfaces between domains comprising different formulations. However, Discontinuous Galerkin approximations lead to larger systems of equations with extended bandwidth and, consequently, demanding greater computational capacity.

The application of H(div)-conforming spaces revealed to be an efficient form of approximating velocity fields in Stokes, Darcy and coupled Stokes-Darcy equations, leading to optimal convergence rates, even when velocity and pressure spaces are defined by polynomials of the same degree, as in the schemes using

- → Q d RT (k) Q d k and - → Q div RT (k) Q d
k space configurations based on quadrilateral meshes. Furthermore, the former one has the additional advantage of being strongly conservative, and to allow natural enforcement of Beavers-Joseph-Saffman interface condition in coupled flows.

Figure 1 :

 1 Figure 1: Stokes problem: influence of the penalization parameter γ| E = γ0k 2|E| on the accuracy of u using symmetric (S) and non-symmetric (N-S) methods for space configurations-→ Q d k Q d k-1 ,

Figure 2 :

 2 Figure 2: Stokes problem: errors in u (left side) and in p (right side), measured in L 2 -norms, obtained with formulations based on square meshes with spacing h, for k = 2 (top side) and k = 3 (bottom side): symmetric (S) and non-symmetric (N-S) DG methods of types -→ Q div RT (k) Q d k , and

1 Figure 3 :

 13 Figure 3: Stokes problem: L 2 -errors in u (left side) and p (right side), for H 1 -conforming methods with stable space configurations -→ Q c k Q c k-1 ,

  y) cos x v(y) sin x , p S = sin x sin y, where v(y) = 1 π 2 sin 2 (πy) -2, u D = (e -y -e y ) cos x -(e -y + e y ) sin x , p D = (-e -y + e y ) sin x.

Figure 4 :

 4 Figure 4: Stokes-Darcy problem: convergence history for u (left side), p (right side), on Ω S (top side) and on Ω D (bottom side), using spaces -→ Q div RT (k) Q d k , k = 1, 2 and 3, with symmetric formulation on Ω S .

• Step 1 : 2 [ 8 ]• Step 2 :

 1282 Consider the domain Ω = Ω 1 ∪ Ω 2 (with 8 × 4 size) illustrated in the left side of Fig.5, where a Newtonian fluid (µ = 1) is governed by the Stokes equation. On the right part Ω 2 , the domain is populated with obstacles of radius r = 0.25, where the boundary condition u = 0 is enforced around them (i.e., on ∂Ω f 2 ). At the left border, u • n = ûn = 1, and p n = (2µD( u) -p I) • n = 0 on the right side. Elsewhere, no-penetration condition ûn = 0 is specified. Numerical resolution with H(div, Ω)-conforming formulation is performed, with space configurations of type -→ P div BDF M (3) P d based on the adaptive triangular mesh shown in Fig.6. Coupled Stokes-Darcy problem with domains Ω S = Ω 1 and Ω D = Ω 2 , and same boundary conditions as in Step 1, according to illustration in the right side of Fig.5. The permeability coefficient for Darcy's domain (κ num ) is obtained from the numerical solution of pressure of Step 1, and it is also verified by an homogenization procedure documented in[START_REF] Tamayol | Analytical determination of viscous permeability of fibrous porous media[END_REF]. The problem is approximated by the combined formulation with space configuration-→ V div Ψ d of type -→ Q div RT (1)Q d 1 on both sides. For this simulation, the mesh is much more simpler: exact results are obtained by using a mesh with 4 × 2 quadrilateral elements.

Figure 5 :

 5 Figure 5: Domains and boundary conditions: Stokes flow with obstacles solved in Step 1 (left side), coupled Stokes-Darcy flow solved in Step 2 (right side).

Figure 6 :

 6 Figure 6: Step 1: mesh geometry (1255 nodes and 2540 elements) for the Stokes flow with obstacles.

Figure 7 :

 7 Figure 7: Fluid pressure fields: Stokes flow with obstacles solved in Step 1 (left side), coupled Stokes-Darcy flow solved in Step 2 (right side).

Figure 8 :

 8 Figure 8: Fluid pressure profiles at y = 2 for Step 1 (left side), and Step 2 (right side).

Table 2 :

 2

	8.82E -5	4.45E -4	8.79E -5	4.43E -4

Stokes problem: errors in u and p, measured in L 2 -norms, obtained with formulations based on square meshes with spacing h = 1/N , and using k = 2 (left side) and 3 (right side): symmetric (S) and nonsymmetric (N-S) DG method of types -

Table 3 :

 3 Stokes problem: L 2 -errors in ∇• u for the H 1 -conforming method of types

		3.70E -1	5.14E -2	4	4.38E -14	4.36E -14
	8	9.97E -2	6.78E -3	8	7.06E -14	6.62E -14
	16	2.54E -2	8.93E -4	16	9.91E -14	9.91E -14
	32	6.37E -3	1.28E -4	32	1.66E -13	1.64E -13
	64	1.60E -3	2.25E -5	64	3.33E -13	3.20E -13

Table 4 ,

 4 κ est and κ num values are compared for different meshes with obstacles of radius r.

	r	0.25	0.30	0.35	0.40
	κest	0.0250 0.0130 0.0050	0.0020
	κnum	0.0212 0.0115 0.0054	0.0019
	∆pnum	188.44 348.57 738.43 2067.32

Table 4 :

 4 Comparative study of permeability coefficients κ num (corresponding to ∆p num ) and κ est

Acknowledgments

The authors thankfully acknowledge financial support from CAPES -Brazilian Federal Agency for Post-graduate Education -(grant 01P-4376/2015), CNPq -the Brazilian Research Council (grants 305425/2013-7, and 304029/2013-0), and to FAPESP -Research Foundation of the State of São Paulo, Brazil (grant 2016/05155, and 2016/13637-5).

Weak formulations for the Stokes-Darcy problem

The main goal here is to show how the well known capacity of divergence-conforming spaces for flux representation in mixed methods for Darcy's flows can be combined with the different schemes for Stokes flows illustrated in the previous section to design different stable and accurate weak formulations for the coupled Stokes-Darcy problem.

Let T S and T D be shape-regular partitions of Ω S and Ω D , respectively, and assume they match along the interface Γ SD . This means that T = T S ∪ T D form a partition of Ω S ∪ Γ SD ∪ Ω S . Given pressure and velocity approximation spaces Ψ and -→ V based on the partition T , the mixed formulation searches for a pair of functions

The bi-linear forms for the coupled problem can expressed as

where the components a S ( u S , ϕ S ) and b S ( u S , ϕ S ), coming from the Stokes part, are obtained from the corresponding ones presented in the previous section, with some minor modifications. For instance, the integrals are restricted to the elements Ω e ∈ T S , to the edges in the skeleton Γ S that are not included in the interface Γ SD , or to the internal edges in Γ int,S . Following the unified presentation of the previous section for different schemes to approximate Stokes flows, three kinds of approximation spaces

V div can be considered for the approximate velocity fields on coupled Stokes-Darcy problems, depending on the continuity constraint required for the Stokes velocity:

For the velocity space configurations -→ V d,div and -→ V div , it is convenient to take local polynomial spaces -→ V e of the same type in both sides, for Ω e ∈ T S and for Ω e ∈ T S . But they should be well balanced with the corresponding pressure approximation space Ψ = Ψ d by the compatibility property ∇ • -→ V e = Ψ e . The methodology using vector spaces of type -→ V d,div has been introduced and analyzed in [START_REF] Rivière | Locally conservative coupling of Stokes and Darcy flows[END_REF] by using local approximations of type P k P k-1 . Then the method has been extended in