O L Cruz-González

R Cornaggia

S Dartois

R Brenner

Accounting for spatial distribution in mean-field homogenization of particulate composites

Keywords: microstructures, particulate reinforced composites, homogenization, spatial distribution, Voronoï diagram

Several analytical mean-field homogenization methods, which take into account the particle volume fraction, shape and orientation are readily available to estimate the effective properties of particulate composites. Models have also been proposed to account for the spatial distribution of the particles. The classical Ponte-Castañeda and Willis (PCW) model is based on a parametrization of the statistical distribution law, while the Interaction Direct Derivative (IDD) model associates a matrix cell to each inclusion, representative of close interactions. In the literature, the use of the IDD is commonly reduced to the particular case of the classical Mori and Tanaka (MT) scheme or to the aforementioned PCW model.

The present study introduces an original approach to calibrate the IDD model, for 2D linear conductivity, based on representative images of the microstructure. The links between the models and the range of validity of the IDD model are discussed. Besides, an "IDD-based" PCW model and a two-step scheme are proposed for situations where the IDD estimate is inconsistent (lack of major symmetry). Finally, an image analysis method using Voronoï diagrams is implemented to define the cells associated to each inclusion and supply the models. The method is validated by comparisons between the obtained IDD and PCW estimates, the Mori-Tanaka (MT) model and benchmark full-field numerical simulations.

Accounting for the inclusion distribution is seen to lead to better estimates, both qualitatively (by capturing anisotropic behaviors due to the sole distribution) and quantitatively.

Possible extensions to elastic composites are discussed.

Introduction

The ongoing development of more efficient and sustainable composite materials is a current phenomenon across multiple industries. Nowadays, a variety of composites are created with the aim of realizing a spectrum of objectives, such as weight reduction, enhanced permeability, augmented sound or thermal insulation, and amplified ductility or strength [START_REF] Khan | A state-of-the-art review on particulate wood polymer composites: Processing, properties and applications[END_REF][START_REF] Rajoriya | Thermal conductivity estimation of high solid loading particulate composites: A numerical approach[END_REF][START_REF] Dib | Generalized modeling of the effective thermal conductivity of particulate composites[END_REF].

During this process, the effective behavior of the materials is often evaluated through homogenization methods. These methods are used to determine the macroscopic response of the material when it is exposed to external loads such as thermal or mechanical forces at a significantly larger scale, which is well separated from the microstructural scale. This approach enables to investigate the behavior of materials with varying microstructures and gain a comprehensive understanding of their mechanical and thermal properties at the macroscopic level.

In the literature, several micromechanical homogenization methods are available. On the one hand, numerical full-field approaches [START_REF] Bluthé | Contribution of FE and FFT-based methods to the determination of the effective elastic and conduction properties of composite media with flat inclusions and infinite contrast[END_REF][START_REF] Schneider | A review of nonlinear FFT-based computational homogenization methods[END_REF][START_REF] Lucarini | FFT based approaches in micromechanics: fundamentals, methods and applications[END_REF] provide the "exact" effective properties of specific microstructures, up to numerical approximation. As a downside, they may be computationally expensive, especially when studying the influence of microstructural variations by performing multiple evaluations. On the other hand, mean-field models are based on the estimation of average fields of interest as stresses, strains, heat flux, intensity, etc., in the constituents of the composite [START_REF] Zaoui | Continuum micromechanics: Survey[END_REF]. In the context of matrix-inclusion composites, approximating the inclusion shapes by ellipsoids (ellipses in 2D) allows to use the Eshelby method [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] to produce such estimates. Thus, a quick calculation of the effective properties from a representative image or statistical descriptors of the microstructure can be performed.

As stated in [START_REF] Zheng | An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution[END_REF], an ideal mean-field estimation scheme for the effective properties of composites should adequately account for the influence of the inclusion distribution and the interaction between inclusions and their immediately surrounding matrix. In particular, the spatial distribution emerges as a relevant statistical descriptor of a composite, and may depend on manufacturing process parameters such as, injection speed, mixing process, etc., especially for medium and high volume fractions of fibers [START_REF] Nguyen | Fiber length and orientation in long-fiber injection-molded thermoplastics -part I: modeling of microstructure and elastic properties[END_REF][START_REF] Müller | Homogenization of linear elastic properties of short-fiber reinforced composites -A comparison of mean field and voxel-based methods[END_REF]. The Mori-Tanaka (MT) model [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF] is a well-established approach for analyzing multi-phase composites. It incorporates the effects of inclusion interaction and is applicable for a wide range of inclusion geometries, including both isotropic and anisotropic materials. Nevertheless, the MT model does not account for the spatial distribution of inclusions. Conversely, there exist advanced mean-field models which are capable of addressing this limitation of the MT model, such as the Ponte-Castañeda and Willis (PCW) model [START_REF] Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF], which relies on a statistical description of the phase distribution, and micromechanical coated-inclusion approaches [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinders models[END_REF][START_REF] Hori | Double-inclusion model and overall moduli of multiphase composites[END_REF][START_REF] Aboutajeddine | The double-inclusion model: a new formulation and new estimates[END_REF] that lead to the derivation of the Interaction Direct Derivative (IDD) model [START_REF] Zheng | An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution[END_REF][START_REF] Du | A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution[END_REF].

Based on the above, this study introduces a novel approach for calibrating and utilizing the IDD model and a proposed IDD-based PCW model in the homogenization of two-dimensional (2D) linear conductive composite materials. The composites consist of a matrix reinforced by elliptical inclusions, both isotropic materials. Four different modeling approaches are applied: MT, PCW, IDD, and a FFT-based full field approach (FF). While MT and FF are widely used in the literature with established parameter selection methods, the situation is somewhat distinct for the classical PCW, which is also frequently employed, and for IDD. The selection of cell parameters from a given microstructure remains an open question, despite their use in previous studies [START_REF] Charpin | Estimating the poroelastic properties of cracked materials[END_REF][START_REF] Müller | Homogenization of linear elastic properties of short-fiber reinforced composites -A comparison of mean field and voxel-based methods[END_REF][START_REF] Hessman | On mean field homogenization schemes for short fiber reinforced composites: Unified formulation, application and benchmark[END_REF][START_REF] Suarez-Afanador | Effective thermo-viscoelastic behavior of short fiber reinforced thermo-rheologically simple polymers: An application to high temperature fiber reinforced additive manufacturing[END_REF][START_REF] Zhou | An improved mean-field homogenization model for the three-dimensional elastic properties of masonry[END_REF].

The aim of this work is to accurately estimate the effective behavior of 2D conductive composite materials by leveraging the IDD method and an IDD-based PCW approach, utilizing microstructure images. Moreover, the research strives to determine under which conditions these methods offer better outcomes compared to MT. FF simulations results will be used as references. It is recalled that the IDD estimate is inconsistent in particular situations (lack of major symmetry). In these specific cases, an IDD-based PCW model and the two-step homogenization approach from [START_REF] Pierard | Mean-field homogenization of multi-phase thermoelastic composites: a general framework and its validation[END_REF] are proposed to overcome this issue. The main novelty of the paper lies in proposing an image processing technique using Voronoï diagrams that permits the selection of cells associated to each inclusion in the IDD scheme. Furthermore, it is proposed to define the identical distributional cells and the corresponding Hill polarization tensor in the PCW scheme in such a way that this method connects directly to the IDD approach and includes averaged information coming from the Voronoi diagrams.

Finally, we present our findings for a variety of artificial microstructures chosen to exhibit "distributional" effects and to underline the relevance of the presented models.

The manuscript is organized as follows. In Section 2, we give an overview of the homogenization methods in the context of the conductivity problem, where the existing links between MT, PCW and IDD models are discussed. In Section 3, we present our main proposal concerning the properties and practical implementation of the IDD model, and notably the determination of the cells representative of the spatial distribution. Section 4 describes the methods, modeling methodology and numerical tools used for the purposes of the manuscript. Then, in Section 5, we explore 3 different microstructural frameworks that hold fundamental significance, and perform comparisons between the homogenization schemes. Finally, in Section 6 we summarize our findings and discuss the limitations and possible extensions of our work. Well-known useful formula and auxiliary results are gathered in appendices for completeness.

Homogenization of conductive composites

Let us consider a heterogeneous, linear conductive material occupying a domain B (see Figure 1 (a)). Analogous to the measurement of the macroscopic properties of a material from a representative test specimen, here, it is considered a representative volume element (RVE) V (see Figure 1 (b)) on the microscopic level, which has to be representative of the entire material.

Focusing on matrix-inclusions composites, V M = V 0 refers to the matrix domain and

V I = ∪ N α=1 V α represent the inclusions domain with N ∈ N.
The interfaces between inclusions and matrix, denoted I , are assumed to be perfect.

In the absence of heat sources, the temperature field u and the heat flux q satisfy the thermal equilibrium equation and transmission conditions across interfaces: where n denotes the normal vectors to the interface I , and • denotes the jumps across the interface I .

∇ • q = 0 in V \ I , u = 0 and q • n = 0 on I , (1)

100

Additionally, the temperature u(x) is linked to the heat flux q(x) by Fourier's law,

q(x) = -k(x)g(x), (2)
where g(x) := ∇u(x) represents the local intensity field, and k(x) refers to the second-order conductivity tensor which is symmetric and positive definite.

The macro-heat flux q and macro-intensity g which characterize the state of the material at the macroscopic level are defined as the averages of their microscopic counterparts in the RVE, that is,

q := 1 |V | V q(x)dV , g := 1 |V | V g(x)dV . (3)
The aim of the homogenization is to obtain the effective behavior that describes the composite at the macro scale based on the information at the micro scale, i.e. the constant effective conductivity tensor k eff that relates the macroscopic flux and intensity as follows:

q = -k eff g . (4)
This effective tensor can be computed numerically by solving elementary problems on the RVE to determine the local fields (hence the name full-field for this approach) for various representative load cases, as precised for completeness in Appendix A. These full-field problems can be solved using several numerical methods, e.g. finite elements [START_REF] Hussein | Micromechanics based FEM study on the mechanical properties and damage of epoxy reinforced with graphene based nanoplatelets[END_REF][START_REF] Bluthé | Contribution of FE and FFT-based methods to the determination of the effective elastic and conduction properties of composite media with flat inclusions and infinite contrast[END_REF] or FFT-based methods [START_REF] Schneider | A review of nonlinear FFT-based computational homogenization methods[END_REF][START_REF] Lucarini | FFT based approaches in micromechanics: fundamentals, methods and applications[END_REF]. However (i) these computations can be costly (not in the 2D conductivity examples provided in the present paper but for instance for 3D elasticity and large RVE, or when many microstructures should be evaluated as in the case of parametric sensitivity studies) and (ii) require a perfect knowledge of the microstructure in the RVE (i.e. the conductivity field k(x)), while one often only has access to statistical information about a composite.

In contrast with these full-field methods, mean-field methods provide only estimates of the effective properties, by approximating the means of the fields in each constitutive phase for piecewise-homogeneous composites. These semi-analytical estimates are much faster to compute and thus emerge as an alternative to FF computations. They are presented now.

Mean-field homogenization

Let us consider now a multi-phase composite with piecewise constant properties, so that the conductivity field can be stated as follows:

k(x) = N α=0 k α χ α (x), (5)
where k α is the conductivity tensor of the α-th phase, and χ α refers to the characteristic function of the domain V α .

The aim of the mean-field methods is to obtain estimates of the mean values q α and g α in the different phases, given the overall means q and g . Here, the notation • α refers to the mean over the α-th phase, namely

• α := 1 |V α | Vα (•)dV α , (6)
so that • = α=0 c α • α and c α is the volume fraction of the phase V α .

Since the thermal equilibrium problem is linear, there exist the so-called second order localization tensor A and the mean tensors A α = A α that link microscopic and macroscopic quantities:

g(x) = A(x) g , and

g α = A α g . (7)
The following property is readily fulfilled:

A = N α=0 c α A α = I. (8)
At this point, if the mean localization tensors are known, one can compute the average flux as:

q = N α=0 c α q α = N α=0 c α (-k α g α ) = - N α=0 c α k α A α g , (9)
and by analogy with Eq. ( 4), the effective conductivity is given by:

k eff = N α=0 c α k α A α = kA . (10)
Additionally, for matrix-inclusions composites, Eq. ( 10) is commonly rewritten in such way that it only depends on the inclusion localization tensors using the property [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] to remove the mean matrix localization tensor A 0 :

k eff = k 0 + N α=1 c α (k α -k 0 )A α . (11)
Finding the exact values of mean localization tensors A α would require to solve local equilibrium problems such as (A.1a)-(A.1d). Instead, mean-fields methods look for approximations of these tensors, using hypotheses on the phase geometries and spatial distributions.

Considered mean-field schemes

Eshelby's equivalent inclusion method was proposed by Eshelby in [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] to determine the elastic solutions of a single ellipsoidal inclusion embedded in an infinite matrix with uniform external loading. This micromechanical estimation has been the backbone of most meanfield methods. In particular, when the volume fraction of the inclusions in a composite is very low, the interactions between the inclusions can be considered negligible and the dilute approximation arises naturally from Eshelby's theory, using the formula [START_REF] Müller | Homogenization of linear elastic properties of short-fiber reinforced composites -A comparison of mean field and voxel-based methods[END_REF] with the inclusion localization tensor computed as if it was isolated in the matrix, namely

A (incl) α = I + P (incl) α (k α -k 0 ) -1 , (12)
where P (incl) α stands for the Hill tensor of the α-th inclusion that accounts for its geometry: see [START_REF] Parnell | The Eshelby, Hill, moment and concentration tensors for ellipsoidal inhomogeneities in the Newtonian potential problem and linear elastostatics[END_REF] and Appendix B.2 for specific expressions.

The following is a brief summary of the mean-field methods of interest for the present work. We extensively rely on the formalism chosen by [START_REF] Müller | Homogenization of linear elastic properties of short-fiber reinforced composites -A comparison of mean field and voxel-based methods[END_REF][START_REF] Hessman | On mean field homogenization schemes for short fiber reinforced composites: Unified formulation, application and benchmark[END_REF] to write the effective stiffness tensors of elastic composites. Here, we adapt these expressions for the purposes of the conductivity problem.

Mori-Tanaka approximation

In the Mori-Tanaka scheme [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] as reformulated by Benveniste [START_REF] Benveniste | On the effective thermal conductivity of multiphase composites[END_REF][START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF], the field in the matrix at a sufficient distance from an inclusion is approximated by the constant value of its mean. As a result, the existence of further inclusions is encoded in the mean field of the matrix, and thus the method takes into account the particle interactions. The mean field of the inclusion in this model is linked to the mean field of the matrix by Eshelby's localization tensor [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] instead of the total mean field.

Based on these considerations, the mean localization tensor of the α-th inclusion takes the following form:

A (MT) α = A (incl) α c 0 I + N β=1 c β A (incl) β -1 . (13)

Ponte-Castañeda and Willis estimate

The Ponte Castañeda and Willis scheme was proposed in [START_REF] Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF] to address the effective behavior of particulate composites with ellipsoidal spatial correlation of inclusions positions.

This variational Hashin-Shtrikman type estimate [START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF] accounts for inclusion shape and spatial distribution independently. In the following, we will consider the restriction of the PCW model to the case of identical shape of spatial correlations between inclusions. Readers are referred to [START_REF] Duan | Effective conductivities of heterogeneous media containing multiple inclusions with various spatial distributions[END_REF] for the application of the PCW model to conductive composites.

The localization tensors for this scheme is given as follows,

A (PCW) α = A (incl) α I -P (cell) N β=1 c β (k β -k 0 ) A (incl) β -1 , (14)
where P (cell) stands for the Hill tensor of the identical distributional cells that characterizes the spatial distribution of the inclusions.

Interaction Direct Derivative model

The Interaction Direct Derivative model was introduced by Zheng and Du in [START_REF] Zheng | An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution[END_REF][START_REF] Du | A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution[END_REF] from a three-phase approach: a matrix cell of ellipsoidal geometry surrounds each inclusion and in turn is surrounded by an unbounded medium with effective properties. In [START_REF] Deng | Interaction models for effective thermal and electric conductivities of carbon nanotube composites[END_REF], the IDD model applied to conductive composites is presented.

By following this considerations, the localization tensors

A (IDD) α
for the IDD model takes the form

A (IDD) α = I - N β=1 c β (k β -k 0 ) A (incl) β P (cell) β -1 A (incl) α , (15)
where

P (cell) β
represents the Hill tensor of the β-th matrix cell.

To summarize, although MT scheme takes into account the inclusion interaction, it fails to account for the inclusion distribution and the interaction between inclusion and their immediate surrounding matrix material. On the other hand, PCW and IDD estimates do not present the aforementioned shortcomings and capture this information through the matrix cells.

Relations between IDD, MT and PCW models

In order to study the connections between the mean-field approaches IDD, MT and PCW, a unified formulation should be considered, as proposed by [START_REF] Hessman | On mean field homogenization schemes for short fiber reinforced composites: Unified formulation, application and benchmark[END_REF], based on the links established earlier between these models in [START_REF] Du | A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution[END_REF].

In this vein, after transformations on the localization tensor A (MT) α of Eq. ( 13), the following result is obtained:

A (MT) α = A (incl) α I - N β=1 c β P (incl) β (k β -k 0 )A (incl) β -1 . (16)
Then, combining the three approximations, ( 14), ( 15) and ( 16) of mean localization tensors with the definition (11) of the effective conductivity, the three estimates are given in Table 1, with a transpose on the IDD estimate to stress the similarity of expressions. Indeed, the IDD estimate may not be symmetric, an important issue addressed in detail in Section 3.1 below.

Approaches

Effective conductivity tensors MT k

(MT) eff = k 0 + N α=1 c α (k α -k 0 ) A (incl) α I - N β=1 c β P (incl) β (k β -k 0 )A (incl) β -1 , PCW k (PCW) eff = k 0 + N α=1 c α (k α -k 0 ) A (incl) α I -P (cell) N β=1 c β (k β -k 0 ) A (incl) β -1
,

IDD T k (IDD) eff = k 0 + N α=1 c α (k α -k 0 ) A (incl) α I - N β=1 c β P (cell) β (k β -k 0 ) A (incl) β -1
.

Table 1: Unified formulation of the effective conductivity properties. The IDD estimate is transposed to stress the similarity of expressions.

As observed, the expressions in Table 1 are very similar, and when the IDD estimate is symmetric, i.e. k

(IDD) eff = T k (IDD) eff
, a link among IDD, MT, and PCW can be provided, as discussed by [START_REF] Du | A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution[END_REF][START_REF] Hessman | On mean field homogenization schemes for short fiber reinforced composites: Unified formulation, application and benchmark[END_REF] and recalled now for completeness. In Figure 2 (a), the capability of IDD to account for inclusions and cells with different aspect ratios and orientations is illustrated. In contrast, the MT approach (see Figure 2 (b)) considers that inclusions and cells respectively share the same aspect ratio and orientation, which is achieved by setting

P (cell) β = P (incl) β
in the IDD scheme. Besides, the PCW scheme is obtained by replacing 1).

P (cell) β = P (cell) (see Table
Among the studied models, the IDD approach offers morphological criteria by utilizing cells that consider the spatial distribution of inclusions and their interactions with the immediately surrounding matrix. Such cell parameters should be chosen or derived from microstructural analyses. Some propositions are done by [START_REF] Zhou | An improved mean-field homogenization model for the three-dimensional elastic properties of masonry[END_REF][START_REF] Suarez-Afanador | Effective thermo-viscoelastic behavior of short fiber reinforced thermo-rheologically simple polymers: An application to high temperature fiber reinforced additive manufacturing[END_REF], and a systematic approach is presented below.

Remark 1. The coated-inclusion approach proposed prior to IDD by [START_REF] Hori | Double-inclusion model and overall moduli of multiphase composites[END_REF] results in very similar estimates than IDD, and also reduces to the MT and PCW models in special cases [START_REF] Aboutajeddine | The double-inclusion model: a new formulation and new estimates[END_REF]. The question whether IDD and this coated-inclusion model are identical in all cases has not been answered yet, at the best of our knowledge, but is out of the scope of the present paper.

Practical use of the IDD model

In this section, we address two key issues that have to be tackled to use effectively the IDD model, that is (i) how to account for its (possible) lack of symmetry and (ii) how to choose the matrix cells corresponding to each inclusion.

Addressing symmetry issues

In relation to the MT scheme, several studies have addressed the problem of loss of symmetry appearing in the effective tensor. For example, [START_REF] Ferrari | Asymmetry and the high concentration limit of the Mori-Tanaka effective medium theory[END_REF] dealt with two-phase composites and showed that if the composite is reinforced with isotropic inclusions of any morphology or with perfectly aligned fibers of any material symmetry, these are sufficient conditions to guarantee the symmetry of the MT tensor. Note that the term "phase" stands for inclusions of assigned geometry and constituent material, regardless of orientation. In addition, the authors in [START_REF] Benveniste | On diagonal and elastic symmetry of the approximate effective stiffness tensor of heterogeneous media[END_REF] focused their study on multi-phase composites, and proved that the MT estimates is symmetric only for those multi-phase composites where all phases have similar shape and the same orientation. By contrast, it is stressed out that the PCW estimate always has the correct symmetry for the effective conductivity tensor [START_REF] Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF][START_REF] Duan | Effective conductivities of heterogeneous media containing multiple inclusions with various spatial distributions[END_REF].

Finally, regarding IDD, we proceed in a similar manner as was discussed in [START_REF] Hessman | On mean field homogenization schemes for short fiber reinforced composites: Unified formulation, application and benchmark[END_REF] and rewrite the IDD expression given in Table 1 as follows:

T k (IDD) eff = k 0 +   N α=1 c α (k α -k 0 ) A (incl) α -1 -P (IDD)   -1
with:

P (IDD) = N β=1 c β P (cell) β (k β -k 0 ) A (incl) β N α=1 c α (k α -k 0 ) A (incl) α -1 , (17)
Note that this is almost exactly the PCW expression (see Table 1), only with P (cell) replaced 200 by P (IDD) . The symmetry of the IDD model is thus guaranteed when this tensor

P (IDD) is symmetric, since (k α -k 0 ) A (incl) α
is always symmetric. This is in particular the case in the following configurations:

1. the MT model is used (Figure 2 = A (incl) . In this case, one obtains

P (IDD) = 1 c I N β=1 c β P (cell) β , with c I = N α=1 c α (18)
which is symmetric as a sum of symmetric Hill's tensors;

4. all inclusions are circular (in 2D) or spherical (in 3D) and made of isotropic but not necessarily identical materials, with conductivity tensors k α = k α I. In this case, the localization tensors are isotropic (proportional to identity) and the expression [START_REF] Du | A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution[END_REF] of P (IDD) again becomes a weighted sum of symmetric Hill tensors. In 2D, with the expressions given in Appendix B.2 one obtains:

P (IDD) = N β=1 c β P (cell) β 2(k β -k 0 ) γ β + 1 N α=1 c α 2(k α -k 0 ) γ α + 1 -1 (2D case) ( 19
)
where γ α = k α /k 0 are the conductivity contrasts. This last case is of interest when studying the behavior of composites reinforced with two (or more) kinds of long, unidirectional fibers by reducing the analysis in the (2D) transverse plane e.g. [22, Sect.

5.2.5]

; or 3D composites containing both seemingly-spherical reinforcements and voids.

The third and fourth cases above are configurations where IDD can be used "as it is", without reducing to one of the two other models. The third one will be used in the ensuing numerical examples for comparisons. To address the cases where symmetry is not guaranteed, while still accounting for the relevant information given by the individual cells, we now propose two methods.

Proposed IDD-based PCW model

A primary consideration is to make use of the PCW model. The question then is the choice of the distributional cell and the corresponding tensor P (cell) . Based on the third case and on expression [START_REF] Charpin | Estimating the poroelastic properties of cracked materials[END_REF] just above, we propose the following choice for this tensor:

P (cell) := 1 c I N α=1 c α P (cell) α . (20)
In this way, IDD and PCW models coincide in cases 2 and 3 discussed above (identical cells and identical inclusions), and the symmetric estimate whose distributional tensor is given by ( 20) can be used in other configurations. This distributional tensor is used for the numerical simulations of the PCW model in Section 5 below.

Two-step homogenization process

Another approach for dealing with the symmetry issue of the MT model was proposed in [START_REF] Pierard | Mean-field homogenization of multi-phase thermoelastic composites: a general framework and its validation[END_REF] and followed by many afterwards, e.g. [START_REF] Tian | Multi-scale and multi-step modeling of thermal conductivities of 3D braided composites[END_REF]. A two-step homogenization scheme is applied to compute the effective properties combining MT and the upper and lower bounds of Voigt and Reuss, respectively, namely

k (Voigt) eff := N α=0 c α k α , (21a) k
(Reuss) eff := N α=0 c α (k α ) -1 -1 . (21b)
Here, instead of using MT, we adopt the IDD scheme for the first homogenization step.

Figure 3 describes the two-step homogenization process. The domain is first decomposed into as many regions as there are families of identical inclusions. These regions behave as a two-phase material where the volume fraction of inclusions in each grain is the same as the total volume fraction of inclusions in the original composite. As a result, the symmetry of the IDD tensor is guaranteed in all the regions, as discussed in Section 3.1 above. Then, in the first homogenization step we apply the IDD scheme, obtaining the effective properties of each grain individually. The last step is to homogenize over all the grains by means of Voigt and Reuss bounds. We will refer these approaches as IDD-Voigt and IDD-Reuss, respectively.

Remark 2. Other approaches could be used for the second step, as long as they preserve the symmetry of the estimate, e.g. the PCW model, or a self-consistent model. Considering that the geometric features of the microstructure are accounted for by the IDD step, we chose the Voigt and Reuss estimates for the second one for simplicity.

Cells selection via inertia equivalence and Voronoï diagrams

Until now, a comprehensive discussion of the IDD scheme has been provided. However, beyond the particular case where IDD is reduced to MT (see Figure 2), the lack of clarity on how to select the matrix cells (or cell, for the IDD-based PCW model) still persists. In this regard, we propose a methodology to derive the required information of the cells from a microstructural analysis of a representative image. This data will serve not only as input for IDD, but also for applying the IDD-based PCW scheme as proposed in Section 3.1.1 above.

The approach represents one of the main novelties of the present work.

More precisely, we look at the practical implementation of the inertia equivalence procedure proposed by [START_REF] Du | A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution[END_REF] (see Figure 4) but applied only for inclusions distributed on periodic lattices [START_REF] Charpin | Estimating the poroelastic properties of cracked materials[END_REF], to the best of our knowledge. The idea is to consider regular polyhedrons (polygons in 2D) to partition the matrix domain such that each of these subdomains or cells are made up of an inclusion as a nucleus and an immediate surrounding matrix material region as the atmosphere. Then, since non-ellipsoidal cells are obtained to assemble the space, the procedure follows by assigning a matching ellipsoid to each polyhedral cell. In the procedure, the polyhedron is further substituted by an equivalent ellipsoid that has the same moment of inertia as the region. For instance, if the tessellation polygons are rectangles (e.g. in the case of a bi-periodic distribution of inclusions), then elliptical cells are proportional to the maximal inscribed ellipses.

To determine these polyhedral cells, we quote [START_REF] Du | A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution[END_REF] (with our notation): "For any given inclusion V α , to define the inclusion-matrix cell V

(cell) α

we can introduce a plane between V α and its any neighboring inclusion V β that is perpendicular to the line linking the centers of V α and V β and partitions this line in accordance with the ratio of the radii of V α and V β ".

This description naturally relates to Voronoï diagrams [START_REF] Torquato | Random Heterogeneous Materials[END_REF][START_REF] Pollmann | Modeling and computational homogenization of chloride diffusion in three-phase meso-scale concrete[END_REF][START_REF] Moulinec | A simple and fast algorithm for computing discrete Voronoi, Johnson-Mehl or Laguerre diagrams of points[END_REF]. The Voronoï diagram of a given set of discrete points in a Euclidean space, also called "seeds", is a partition of the space into regions, which contain the points of the space closer to a given seed than to any other. Let X be a metric space with distance function d, then, the Voronoï cell R k , associated with the seed P k is defined as follows,

R k = { x ∈ X | d(x, P k ) ≤ d(x, P j ) for all j = k}. (22)
A natural extension to a multiplicatively weighted Voronoï diagram is given by the following definition of a cell:

R k = { x ∈ X | d(x, P k ) / r k ≤ d(x, P j ) / r j for all j = k}, (23)
where r k and r j are the weights associated with the seeds P k and P j , respectively.

The definition of the multiplicatively weighted Voronoï diagram has a strong analogy with the aforementioned approach given in [START_REF] Du | A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution[END_REF], with the Voronoï weights determined by the inclusions radii. However, as a drawback of using this tessellation, Voronoï cells in a multiplicative scheme exhibit circular arcs instead of straight lines, may be disconnected and 265 have holes. We refer to [START_REF] Bock | Generalized Voronoi Tessellation as a Model of Two-dimensional Cell Tissue Dynamics[END_REF] for a detailed study of the multiplicative case with applications in two-dimensional cellular tissues.

Based on the above considerations, in this paper we propose to use the classical Voronoï diagram as the tessellation method, and thus obtain the polygonal cells that will then be replaced by the equivalent ellipses with the same moment of inertia as the region. Although the it might be worth to explore it as a first approximation. Its performances and drawbacks will be discussed in Section 6.

Methods, modeling and numerical tools

The aim of this section is to discuss the main set of tools and methods that have been used

to produce the numerical comparisons presented in the next section. These tools were implemented in an home-made Python library "PyHom: a Python library for homogenization" whose main functions are represented in Figure 5.

This computational and homogenization-based approach is conceived for calculating the effective conductivity properties of composite materials with complex micro-structures in a 2D framework. The backbone of the library relies on the theory described above, and uses not only existing image processing modules, but also functionality and packages built by the authors. The main interest is using PyHom as a toolbox to study the core issues in the application of IDD and the IDD-based PCW.

Hereafter, we describe in more detail the choices we made for the key step of the process, namely (i) RVE generation, (ii) image analysis and Voronoï diagram generation and (iii) homogenized estimates computation.

Generating representative elements

To ease the comparisons between models, we generate artificial RVE with controlled distributions of inclusions. The two sections below present the methods chosen to generate these RVE.

RSA algorithm to place elliptic inclusions with prescribed distribution

We start by specifying the Random Sequential Adsorption (RSA) algorithm [START_REF] Torquato | Random Heterogeneous Materials[END_REF] more closely and see how it can be applied to build a heterogeneous microstructure in a 2D framework. Let us first recall that in the RSA algorithm, the position of a test particle is generated following a uniform distribution. If the test particle intersects with any of the previously added particles, it is discarded; otherwise, the test particle is added to the collection. Table 2 and Figure 6 present an example of RVE construction. A number N (incl) of elliptical inclusions, whose aspect ratio e (incl) and orientation θ (incl) are specified, should be placed to reach a volume fraction c (incl) (which determines the size of these inclusions).

PyHom

Start END

To obtain a controlled distribution not only due to the inclusion shape, elliptical exclusion zones, which we refer to as RSA cells, are placed first. Aspect ratio e (cell-RSA) and orientation θ (cell-RSA) of these cells are also specified a priori, as well as a security factor to prevent contact between cells, see Figure 6(a). Their size is determined from a size factor that is the ratio between inclusion and cells major axis.

pixels N (incl) c (incl) e (incl) θ (incl) e (cell-RSA) θ (cell-RSA) security factor size factor 256 8 0.05 0.5 -π/4 0.7 π/4 1.3 0.5 As seen on Figure 6(b), the resulting RVE is periodized to facilitate the full-fields computations, i.e. inclusions that cross the edges are duplicated. This also ensures that the chosen volume fraction is reached within the RVE, up to discretization. Finally, an image is generated with a prescribed discretization (256×256 pixels in all the upcoming examples), see Figure 6(c). This image will be the input of both full-fields and image-based mean-fields homogenization methods, as indicated in Figure 5.

Advani-Tucker law for orientation distributions

To generate distribution of non-aligned inclusions, we use the Advani-Tucker law for orientation distributions. In [START_REF] Advani | The use of tensors to describe and predict fiber orientation in short fiber composites[END_REF], this distribution function is proposed to describe an aligned and axisymmetric orientation state in a 3D configuration, see for instance [START_REF] Advani | The use of tensors to describe and predict fiber orientation in short fiber composites[END_REF][START_REF] Suarez-Afanador | Effective thermo-viscoelastic behavior of short fiber reinforced thermo-rheologically simple polymers: An application to high temperature fiber reinforced additive manufacturing[END_REF][START_REF] Bauer | On the dependence of orientation averaging mean field homogenization on planar fourth-order fiber orientation tensors[END_REF]. The explicit form of the distribution function, adapted for a 2D microstructure, is rewritten as follows,

f (θ) = K sin m AT (θ), ( 24
)
where K is a normalization constant and, θ ∈ [0, π] is the inclusion orientation with respect to the chosen privileged direction. This law is characterized by a single parameter m AT . Two extreme cases arise when m AT = 0 for which the law reduces to the uniform distribution, and when m AT → ∞, for which aligned inclusions are recovered.

Derivation and analysis of Voronoï diagram with Python libraries

To build and exploit the Voronoï diagrams needed for distribution analysis, as proposed in Section 3.2, we use two Python libraries: scikit.image [START_REF]Scikit Image -image processing[END_REF] for image analysis, and scipy [START_REF]Scipy -fundamental algorithms for scientific computing in Python[END_REF] for Voronoï diagram generation.

Figure 7 presents the process applied to the example image presented in Figure 6. First, scikit.image [START_REF]Scikit Image -image processing[END_REF] is used to partition the image between subdomains (matrix and inclusions), to extract the characteristics of the inclusions (number, size, aspect ratio, orientation) and to position their centers, see Figure 7(a). Then, these centers are cloned on a 3 × 3 grid, see To analyze this diagram, the same image processing tool can be used: again the library scikit.image [START_REF]Scikit Image -image processing[END_REF] is used to partition the image Figure 7(c) into distinct cells, whose properties are extracted (function regionprops in modulus mesure). These properties notably include the lengths of the minor and major axis of the ellipse with the same moment of inertia as the region, as required by the proposal of Section 3.2. Figure 8 finally shows the elliptic cells resulting from this analysis, along with the inclusions determined in the first step.

FFT-based solver

Reference effective properties are computed using the full-field approach described in Section Appendix A. To solve the elementary problem, a FFT-based solver was implemented in Python. In this work addressing only 2D linear conductivity with moderate material contrasts, the "basic scheme" proposed in the seminal papers of Moulinec and Suquet [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] was chosen for its simplicity of implementation, although many improved FFT-based solvers have been developed in the last decade, see e.g. [START_REF] Schneider | A review of nonlinear FFT-based computational homogenization methods[END_REF][START_REF] Lucarini | FFT based approaches in micromechanics: fundamentals, methods and applications[END_REF].

Numerical results and comparisons

In this section, we analyze different microstructures, and proceed to compute and compare their effective conductivity properties. For the sake of simplicity and focus on spatial distribution, the heterogeneous material B is assumed to be constituted by a matrix reinforced by elliptical inclusions with different aspect ratios and orientations, both the matrix and inclusions being isotropic phases.

The results are presented in terms of the principal components of the conductivity tensors (i.e. their eigenvalues), denoted k eff-I and k eff-II , with the convention k eff-I ≥ k eff-II , corresponding to the conductivities in the main directions (i.e. normalised eigenvectors) n 1 and n 2 . The angle θ eff := (e 1 , n 1 ), called the orientation of the conductivity tensor afterwards, is used to specify these main directions. The notations are represented in Figure 9, where notations for main axes a 1 and a 2 , and orientation θ (incl) of inclusions are also depicted. , respectively.

Case 1: Anisotropic distribution, identical inclusions

The first case of interest concerns identical isotropic inclusions distributed in an anisotropic way in an isotropic matrix. In particular, for circular inclusions, the expected anisotropy of effective conductivity is due to this distribution only. To do so, we generate the microstructure using elliptic RSA-cells by following the process described in Section 4.1.1. The and blue color, respectively. Notice that the Voronoï diagram has resulted in elliptical cells with several shapes and orientations, therefore, in the particular case of IDD model, we must consider as many Hill tensors as available cells (see Eq. (B.1c)).

In addition, Figure 11 (b) shows the scatter plot between aspect ratios e and orientations θ corresponding to the inclusions and cells, using dots and the previously assigned color. We also display the mean value of the properties of the cells by means of a black star, and the pair (θ (cell-PCW) , e (cell-PCW) ) corresponding to the orientation and aspect ratio of the global cell using in the IDD-based PCW scheme, obtained from the definition proposed in Eq. ( 20), and indicated on the graph with a purple square. Following this idea, Figure 12 shows the normalized effective conductivity properties for the microstructure of Figure 10 (b), but varying the aspect ratio e (incl) and orientation θ (incl) of the inclusions, and the conductivity contrast γ. The comparisons are made between MT, IDD and the FF numerical approach: in this case IDD and PCW models coincide using the definition of Eq. ( 20) for the PCW "distribution" tensor.

In Figure 12 (a) we highlight two situations relevant to the analysis. The first one, denoted by (I), corresponds to the value e (incl) = 0.6 and θ (incl) = 0 • , where the MT scheme offers better results than IDD in comparison to the FF numerical approach. Nevertheless, in Figure 12 (b) it is observed that by varying the orientation and fixing the remaining parameters, IDD arises with better agreements. The second one (II) refers to the phenomenon of "geometrical" anisotropy that we mentioned at the beginning of the Section for circular inclusions.

For these circular inclusions, the MT scheme could not capture a possible geometrical anisotropy coming from the distribution of inclusions. In contrast, IDD does identify such anisotropy. The results are shown in Figure 12 (c) for different conductivity contrast γ.

While quantitative accuracy is missing as both models underestimate k eff-I , this qualitative improvement is a major advantage of the IDD model. This capability of capturing some of the anisotropy resulting from the distribution is also illustrated in Figure 13 for elliptic inclusions.

The orientation of the MT effective tensor is again entirely given by the orientation of inclusions while the IDD model captures a small amount of the orientation difference observed on FF results due to the distribution.

Case 2: Approximately isotropic distribution and aligned elliptic inclusions

In contrast to Section 5.1, here, we aim to arrange identical elliptical isotropic inclusions in an approximately isotropic way within an isotropic matrix, so that the anisotropy comes mainly from the orientation and aspect ratio of the elliptical inclusions. For this purpose, we consider circular RSA-cells and the set of parameters shown in Table 4. Figure 14 illustrates a possible RVE that can be constructed under these considerations.

Again, we verified that the chosen RVE is statistically representative by comparing the properties of 10 RVE generated with the same parameters, see Figure C.24 in appendix. The maximum difference of the larger conductivity k eff-I with its mean value over the collection of RVE is this time larger, about 2.5%, but we still consider an individual RVE as representative. The normalized effective conductivity properties of the microstructures obtained by varying the aspect ratio e (incl) and orientation θ (incl) of the inclusions are given in Figure 16. Again, for identical inclusions the IDD and PCW models coincide. As observed in Figure 16 (a), the IDD scheme has a better agreement than MT with the FF numerical approach. Furthermore, in Figure 16 For completeness, we show in Figure 17 the comparison between the orientations of 435 effective tensors for the different schemes, which are clearly given by the orientations of inclusions, as expected for this nearly-isotropic distribution.

Case 3: Non-aligned inclusions

In this section, we proceed to investigate a scenario in which the inclusions are not aligned, i.e. cases where the effective conductivity tensor of the IDD model loses its symmetry. To this end, we consider microstructures in which all inclusions are identical except for their orientation. Two examples are considered: (i) manually chosen orientations and (ii) orientations following the Advani-Tucker law. In both cases, the RSA-cells are chosen identical to each inclusion (slightly magnified by a security factor), and the distribution is therefore only constrained by the inclusion shapes, contrarily to the previous cases.

Inclusion orientation chosen in a given set

For this part, the parameters utilized in the RSA algorithm are presented in Table 5, where θ sample is a set of 30 orientation values manually selected from the set of values {-60, -45, -30, 0, 10, 25, 50, 75} degrees. Figure 18 presents the obtained RVE.

Then, Figure 19 (a) presents the Voronoï diagram corresponding to the inclusions and the inertia-equivalent cells. Figure 19 (b) displays the aspect ratios and orientations for both the inclusions and cells, as well as the mean values for the cell properties and the corresponding pixels N (incl) c (incl) e (incl) θ (incl) e (cell-RSA) θ (cell-RSA) security factor size factor 256 30 0.2 0.5 θ sample e (incl) θ (incl) 1.1 1

Table 5: Elliptic inclusions with different orientations: chosen parameters for RVE generation, as represented in Figure 18.

Inclusion orientation following the Advani-Tucker law

As a final result to illustrate the relevance of applying IDD and IDD-based PCW in the 2D effective conductivity context from a microstructure image, we use the Advani-Tucker (AT) law to generate the orientations of the inclusions as depicted in Section 4.1.2. This law is used to simulate the fibers orientation in short-fibers composites in [START_REF] Nguyen | Fiber length and orientation in long-fiber injection-molded thermoplastics -part I: modeling of microstructure and elastic properties[END_REF][START_REF] Suarez-Afanador | Effective thermo-viscoelastic behavior of short fiber reinforced thermo-rheologically simple polymers: An application to high temperature fiber reinforced additive manufacturing[END_REF] among others. Natural questions are (i) whether and how this orientation distribution may be related to the spatial distribution of fibers (both of them being determined by the manufacturing process in actual composites, and by the RSA process in our numerical study); and (ii) if this spatial distribution is worth considering in mean-field estimates.

In this part, as a preliminary investigation, we study the simplified case of composites reinforced by identical ellipses, of aspect ratio 0.5 (i.e. much more than actual short fibers), oriented following the AT law. Specifically, in each iteration of the calculation, we fix the parameter m AT , generate the orientations of the inclusions, and pass this information to the RSA algorithm to generate 10 microstructures that are used to calculate the normalized effective conductivity properties. Figure 21 presents four of these RVE for increasing values of m AT (and thus increasing alignment of inclusions). more complex anisotropy issues arising from the higher-order tensor formalism and see if the present conclusions are still applicable. Convenient algorithms to generate Voronoï diagrams are also available in 3D e.g. [START_REF] Moulinec | A simple and fast algorithm for computing discrete Voronoi, Johnson-Mehl or Laguerre diagrams of points[END_REF].

2. Classical Voronoï diagrams using only the centers of inclusions were used here for simplicity of implementation. In particular, it enabled the use of the same diagram for microstructures made of "pinned" inclusions varying in orientation, volume fraction, and aspect ratio; and therefore to speed up the computation of the effective properties when generating several artificial test cases. To use additively or multiplicatively weighted Voronoï diagrams as discussed in Section 4.2, it would be necessary to determine the properties of the cells for each microstructure. This step may be in particular necessary to address larger aspect ratio of particles, for which classical Voronoï diagram would lead to cells borders intersecting the inclusions.

3. Of course, an important step towards applications to actual composite would be the use of real microstructure images, either 2D for aligned long-fibers or 3D.

4. Finally, a long-term goal is to rely on relevant statistical descriptors of the spatial distribution e.g. correlation functions [START_REF] Torquato | Random Heterogeneous Materials[END_REF] to parametrize these models, rather than image analysis. These descriptors are theoretically linked with the models already, e.g. in the original definition of the PCW model using the two-point correlation functions of the inclusion phases [START_REF] Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF], but were never extracted from an actual microstructure, to the best of our knowledge. An intermediate step would be to generate artificial images given such descriptors, and then apply the presented methodology on these images.

CRediT authorship contribution statement OLCG: Methodology, Software, Visualization, Formal analysis, Writing -original draft. which is proven to be symmetric using reciprocity identities between the uncoupled problems satisfied by the components U j , that lead to k∇U 12 = k∇U 21 .

Appendix B. Specific for isotropic phases Here we provide the expression of Hill and localization tensors, in 2D, for elliptic inclusions in an isotropic matrix of conductivity k 0 , as given by e.g. [START_REF] Parnell | The Eshelby, Hill, moment and concentration tensors for ellipsoidal inhomogeneities in the Newtonian potential problem and linear elastostatics[END_REF]. Noting (a 1 , a 2 ) the principal directions of the ellipse, see Figure 13 (a), a 1 ≥ a 2 the associated semi-axes and e = a 2 a 1 ≤ 1 its aspect ratio, the Hill tensor is given by:

P (ellipse) =   e k 0 (1+e) 0 0 1 k 0 (1+e)   (a 1 ,a 2 ) . (B.2)
Then, from the definition [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF], when the constitutive material of the inclusion is also isotropic with conductivity k α , the dilute localization tensor is:

A (incl)

Figure 1 :

 1 Figure 1: (a) Macroscale, (b) Microscale: representative volume element (RVE).

Figure 2 :

 2 Figure 2: Diagrams of the unified formulation and the connections between the models.

  (b)) and one of the sufficient conditions given above for the symmetry of the MT estimate is satisfied; 205 2. the PCW model is used (Figure 2(c)); 3. all inclusions are identical except for their size, i.e. have same aspect ratios, orientation and constitutive material, and therefore localization tensors: A (incl) α

Figure 3 :

 3 Figure 3: Two-step homogenization procedure.

Figure 4 :

 4 Figure 4: Inertial equivalence in 2D. a) Partition in polygonal cells. b) A polygonal cell. c) Equivalent elliptical cell with the same geometrical moment. This figure is inspired by Figure 4 of [17].

Figure 5 :

 5 Figure 5: Workflow sketch for the numerical comparison of homogenized estimates for various microstructures.

Figure 6 :

 6 Figure 6: Example of RVE construction using an RSA algorithm. (a) Placement of eight RSA cells: dashed and dotted lines design the cells and their amplification by the security factor, respectively. (b) Placement of inclusions in the central points of the RSA cells and duplication of the inclusions that cross the edges. (c) Final image of the microstructure.

Figure 7 (

 7 Figure 7(b), and the Voronoï diagram is generated by scipy. The regions of interest of this Voronoï diagram will consist of the colored closed polygons, while the outermost regions of the diagram (including unbounded cells) are not involved in the analysis. This generated diagram is then isolated, see Figure 7(c).

Figure 7 :

 7 Figure 7: Voronoï diagram built from image analysis (image built in Figure 6) (a) Center of inclusions. (b) Cloned RVE. (c) Derived Voronoï diagram.

Figure 8 :

 8 Figure 8: Final result of image analysis: inclusions and elliptic cells determined from the Voronoï diagram.

Figure 9 :

 9 Figure 9: Notations for the orientations of inclusions and effective conductivity tensor.

Figure 11 :

 11 Figure 11: Elliptic RSA cells and identical inclusions (circular here): (a) Voronoï diagram, inclusions and respective cells. (b) Scatter plot between aspect ratio and orientation.

Figure 12 :

 12 Figure 12: Elliptic RSA cells and identical inclusions: performance of MT and IDD models in comparison with the FF reference results. (a) Fixed orientation and conductivity contrast, varying aspect ratio of inclusions. (b) Fixed elliptic shape and conductivity contrast and varying orientations. (c) Fixed circular shape, varying conductivity contrast. In all case, the volume fraction of inclusions is 15%.

Figure 13 :Table 4 :

 134 Figure 13: Elliptic RSA cells and identical elliptic inclusions: directions of anisotropy of the effective tensors, as a function of orientations of inclusions.

Figure 14 :

 14 Figure 14: Circular RSA cells and elliptic inclusions: (a) RSA process. (b) 2D representative element of the microstructure.

Figure 15 (

 15 Figure 15 (a) shows the Voronoï diagram corresponding to the microstructure of Figure 14 (b), along with the elliptical inclusions and the different inertially equivalent elliptical cells. Moreover, in Figure 15 (b) we present the scatter plot between aspect ratios e and orientations θ corresponding to the inclusions, cells, the mean value of the cell properties, and the pair (θ (cell-PCW) , e (cell-PCW) ). In contrast to what was observed in Figure 11 (b), while the cells remain elliptical, their aspect ratios approach unity and the orientations exhibit a greater dispersion across the range of values. This is an expected result from the choice of circular RSA-cells to generate the microstructure.

Figure 15 :

 15 Figure 15: Circular RSA cells and elliptic inclusions: (a) Voronoï diagram, inclusions and respective cells. (b) Scatter plot between aspect ratio and orientation.

  (b) we set the value of the largest discrepancy between MT and IDD, i.e. , e (incl) = 0.4, and study the behavior of the properties when varying the orientation of inclusions. The results are still favorable for IDD, proving in this case a clear quantitative gain compared to MT.

Figure 16 :

 16 Figure 16: Circular RSA cells and elliptic inclusions: MT and IDD performance compared to the FF numerical approach for different aspect ratios and orientations of inclusions.

Figure 17 :

 17 Figure 17: Comparison between the directions of anisotropy.

Figure 18 :

 18 Figure 18: Elliptic inclusions with different orientations: (a) RSA process with RSA-cells coinciding with inclusions. (b) RVE of the microstructure. (c) Number of inclusions per chosen orientation.

Figure 19 :

 19 Figure 19: Elliptic inclusions with different orientations: (a) Voronoï diagram, inclusions and respective cells. (b) Scatter plot between aspect ratio and orientation.

Figure 20 :

 20 Figure 20: Elliptic inclusions with different orientations: comparison between the different mean-field approaches and the full-field numerical results. The IDD scheme is assisted by double homogenization (see Section 3.1.2).

Figure 22 presents

 22 Figure 22 presents the mean values and standard deviations of the results obtained using the 10 microstructures for each m AT point. A first observation is that the dispersion among the RVE remains moderate, enabling to clearly discriminate models. For the larger

Figure 21 :

 21 Figure 21: Elliptic inclusions with different orientations following the Advani-Tucker law: examples of RVE for four values of the parameter m AT .

Figure 22 :

 22 Figure 22: Elliptic inclusions with different orientations following the Advani-Tucker law: mean values and standard deviations of the normalized effective conductivity properties obtained using 10 microstructures for each value of the parameter m AT .

RC:

  Funding acquisition, Conceptualization, Methodology, Formal analysis, Writing -original draft. SD: Methodology, Formal analysis, Writing -review & editing. RB: Funding acquisition, Methodology, Formal analysis, Writing -review & editing.

Appendix B. 1 .

 1 Normalized expression of the effective propertiesIf the phases under study are isotropic, we note k α = k α I the conductivity tensor, and the conductivity contrast γ α = k α /k 0 is introduced. The effective properties are normalized as well with respect to the conductivity of the matrix k 0 , and noted keff = k eff /k 0 . From Table1we obtain the normalized estimates:k(MT) eff = I + N α=1 c α (γ α -1) A (incl)of non-dimensional contrasts γ α and normalized Hill tensors P := k 0 P. 560 Appendix B.2. Hill and dilute localization tensors

(a 1

 1 ,a 2 ) , (B.3)

Table 2 :

 2 Parameters for the example RVE represented in Figure6.

classical Voronoï diagram does not seem to conform completely to the authors' original idea,

Acknowledgements

Financial support for this project was provided by the Institute of Materials Science (iMAT) of the Alliance Sorbonne Université through a post-doctoral fellowship to OLCG.

parameters given in Table 3, in the case of circular inclusions, are utilized and result in RVE 365 as the one plotted in Figure 10. pixels N (incl) c (incl) e (incl) θ (incl) e (cell-RSA) θ (cell-RSA) security factor size factor 256 30 0.15 1 0 0.5 0 1.02 0.375 Table 3: Elliptic RSA cells and identical inclusions (circular here): chosen parameters for RVE generation, as represented in Figure 10.

Before going further, the statistical representativity of this RVE is verified by computing Now, the study is focused on a single microstructure, specifically the one shown in Figure 10. In this regard, Figure 11 (a) shows the corresponding Voronoï diagram in the background, along with the circular inclusions and inertially equivalent elliptical cells, highlighted in red

Conclusions and perspectives

This study has addressed the utilization of the Interaction Direct Derivative (IDD) and Ponte-Castañeda and Willis (PCW) mean-field models for estimating the effective conductivity properties of composite materials and accounting for the spatial distribution. We determined configurations for which IDD can be applied "as it is" and other where it lacks the major symmetry. In these later cases, two alternative models, namely an IDD-based PCW model (shading a new light on the links between the two models) and a two-step approach, have been proposed. We then described a methodology that employs Voronoï Diagrams to determine the distributional cells needed to parametrize the IDD model. We developed a home-made Python library that integrates homogenization approaches with image processing tools, enabling the semi-analytical process of simulating the fields of interest.

We applied this framework to determine the effective behavior of 2D conductive composites and compare the performance of IDD and IDD-based PCW against the standard Mori and Tanaka (MT) model and the reference FFT-based full-field numerical simulations.

More precisely, artificial RVE were designed to study the influence of spatial distribution on the effective properties for aligned or non-aligned inclusions. In most of the studied cases, the spatial distribution seems to be an important parameter in the sense that IDD model provides qualitatively or quantitively better results than MT. This in particular striking for non-aligned inclusions, see Section 5.3, where no distribution was imposed a priori, contrarily to the previous cases, and still IDD-inspired models seemed to capture an information not solely contained in the shapes of inclusions (that would be captured by MT).

These are preliminary results on simplified conductive microstructures. This work needs to be further pursued to consider realistic microstructures, representative of e.g. short-fibers composites. At the end, we aim to expand upon the assertion made in [START_REF] Hessman | On mean field homogenization schemes for short fiber reinforced composites: Unified formulation, application and benchmark[END_REF], as stated in the abstract, that "Other models such as the IDD or PCW are of great theoretical importance, but cannot be generally applied for the given material class". To do so, the following points should be explored:

1. The extension to 2D and then 3D elasticity should be performed to explore the much

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Full-field homogenization

Full-field homogenization consists in computing the temperature and flux fields for a basis of elementary macroscopic loads that cover all possible configurations, and then retrieving the effective behavior by computing and relating the means in (3). For conductivity problems, these elementary loads are the unit basis vectors e.g. g = e j if the mean intensity is prescribed. For shorthand notation we regroup them and look for the solution U = {U j } j=1,2

of the following problem:

where I is the second-order identity tensor that regroups the elementary loads: in this The effective conductivity tensor defined by ( 4) is then found to be:

with γ α = k α /k 0 the conductivity contrast. The matrices given in eqs. (B.2) and (B.3) above are the tensors representations in the basis (a 1 , a 2 ) associated with the ellipse, see Figure 9, and classical rotation operators are applied to obtain their counterpart in the fixed basis (e 1 , e 2 ).

565

For circular inclusions, these tensors are isotropic:


Bernard Rolland

Denise Deffains

Patrick Delarue

Stephane Gilles

Antonin Le Campion

None Alain Monnier

Hélène Navier

Alexandre Pichard

Patrice Walczak

Stéphane Gilles

Rémi Perronne

Incidence of four genera of fungi in organic and low-input farming conditions in/on the grain of bread wheat over a 13-year period in France

Keywords: Fungal incidence, Morphological identification, Fusarium, Alternaria, Epicoccum, Microdochium

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

Introduction

The agricultural intensive production model is currently questioned because of agronomic, environmental, public health and sometimes socio-economic issues [START_REF] Van Der Werf | Assessing the impact of pesticides on the environment[END_REF][START_REF] Loyce | Growing winter wheat cultivars under different management intensities in France: A multicriteria assessment based on economic, energetic and environmental indicators[END_REF][START_REF] Bourguet | The hidden and external costs of pesticide use[END_REF]. In this context, reducing the reliance of agriculture on external inputs constitutes a priority to achieve more sustainable models of production.

However, an alternative model such as organic farming has been questioned for several decades about its ability to ensure food safety [START_REF] Magkos | Organic food: buying more safety or just peace of mind? A critical review of the literature[END_REF][START_REF] Lairon | Nutritional quality and safety of organic food. A review[END_REF][START_REF] Garcia | Organic versus conventional food: A comparison regarding food safety[END_REF]. Currently, organic farming is characterized by an increasing consumer demand, which requires paying attention to the sanitary quality of food production, including mycotoxin contamination [START_REF] Brodal | Mycotoxins in organically versus conventionally produced cereal grains and some other crops in temperate regions[END_REF]. This is particularly the case for bread wheat for which grain can be contaminated by several genera of mycotoxin-producing fungi [START_REF] Brodal | Mycotoxins in organically versus conventionally produced cereal grains and some other crops in temperate regions[END_REF], especially by different species of Fusarium causing head blight. This plant disease results from a complex of species belonging to the Fusarium genus [START_REF] Parry | Fusarium ear blight (scab) in small grain cerealsa review[END_REF]. Fusarium head blight (FHB) leads to important economic consequences worldwide, important yield losses and reduction in quality production leading to the destruction of contaminated grain lots, as well as potential contamination by a range of harmful mycotoxins to human and animal health [START_REF] Reddy | An overview of mycotoxin contamination in foods and its implications for human health[END_REF]). These mycotoxins, amongst others deoxynivalenol and zearalenone, are produced by several Fusarium species, especially Fusarium graminearum [START_REF] Reddy | An overview of mycotoxin contamination in foods and its implications for human health[END_REF][START_REF] Dean | The Top 10 fungal pathogens in molecular plant pathology[END_REF].

The mycotoxins are subject to government regulation in the EU, i.e., defined limits in cereals for human consumption (Commission Regulation (EC) No 1881/2006).

Speculation has arisen about the potential higher contamination level of mycotoxins produced by Fusarium species in organic compared to conventional production model, and a set of studies have investigated this issue. However, in the current state of knowledge, organic farming seems usually able to maintain a low level of contamination by mycotoxins, comparable or lower than in conventional farming [START_REF] Brodal | Mycotoxins in organically versus conventionally produced cereal grains and some other crops in temperate regions[END_REF]. The hypothesis that several agricultural practices would be less favourable for the development of Fusarium species causing FHB in organic farming has been proposed. In fact, a low level of nitrogen supply and a crop sequence including fewer cereals reduce the severity of FHB and the level of mycotoxins (Dill-Macky and Jones 2000; [START_REF] Lemmens | The effect of nitrogen fertilization on Fusarium head blight development and deoxynivalenol contamination in wheat[END_REF][START_REF] Heier | Influence of N-fertilization and fungicide strategies on Fusarium head blight severity and mycotoxin content in winter wheat[END_REF][START_REF] Bernhoft | Influence of agronomic and climatic factors on Fusarium infestation and mycotoxin contamination of cereals in Norway[END_REF].

While many studies have focused on the level of mycotoxins in organic and conventional farming (see e.g., [START_REF] Edwards | Fusarium mycotoxin content of UK organic and conventional wheat[END_REF][START_REF] Bernhoft | Influence of agronomic and climatic factors on Fusarium infestation and mycotoxin contamination of cereals in Norway[END_REF][START_REF] Brodal | Mycotoxins in organically versus conventionally produced cereal grains and some other crops in temperate regions[END_REF], only few have considered the incidence of different fungal genera in/on the grain of bread wheat. These latter studies have shown a low incidence of Fusarium species on organically produced grain, while the fungal community was dominated by other fungal genera [START_REF] Lazzaro | Organic vs conventional farming: Differences in infection by mycotoxin-producing fungi on maize and wheat in Northern and Central Italy[END_REF][START_REF] Lenc | Fusarium head blight (FHB) and Fusarium populations in grain of winter wheat grown in different cultivation systems[END_REF].

However, these studies mainly relate to few years and a limited number of varieties, furthermore not identifying the main explanatory factors of the Fusarium incidence in the different production models.

In this study, we monitored, each year over a 13-year period in several sites in France and for a set of bread wheat genotypes, the incidence of four fungal genera on grain of bread wheats from both conventional and organic farming conditions.

More specifically, we addressed the following questions:

(1) What was the incidence of the four genera of fungi on the grain of bread wheat over the period and between sites in both conventional and organic farming systems?

(2) To what extent do year, site, genotype and their interactions explain differences in incidence of the genus Fusarium as well as other fungal genera in/on the grain of bread wheat in both conventional and organic farming systems?

Material and methods

Sample collection

Wheat samples were collected immediately after harvest from experimental field trials from 2007 to 2019 in three areas in France (Brittany, Ile-de-France and Poitou-Charentes) in both conventional and organic farming systems. In our study, conventional farming systems corresponded to low-input systems with no fungicide and growth regulator used, a lower sowing rate and a reduction of nitrogen application (based on N balance sheet minus 60 kgN/ha) but with herbicides used to suppress weeds (see Le Campion et al. 2014 for details).

For each site in low-input farming, when available, a site in organic farming located near the conventional farming site was chosen to limit variation due to large-scale environmental variables such as meteorological conditions (see Online Resource 1 for coordinates). However, no organic farming site was repeatedly available near the low-input site in Poitou-Charentes, so only two sites were monitored in organic farming conditions while three sites were monitored in low-input farming conditions. In each site and each year, four repetitions were carried out for each wheat genotype in plots of 8.5m² using a randomised block design.

Samples were made in conditions of natural contamination on a substantial list of genotypes representative of the winter bread wheat varieties sown in low-input and organic farming over the period considered, including susceptible to tolerant wheat genotypes to FHB. Fusarium graminearum was the main species of Fusarium in France during the period studied. Given the duration of the period studied, the genotypes monitored have been replaced over the years, representing 34 wheat genotypes under low-input conditions and 14 wheat genotypes under organic conditions, for a total of 289 wheat samples collected under low-input and 155 wheat samples collected under organic conditions. Only wheat genotypes monitored for at least two years were retained in this study. The kernels were harvested at maturity with 14% of moisture (0.5-1kg per sample), and sub-samples of 100 kernels were chosen randomly and stored in paper bags at room temperature with 13-15% humidity until analysed for about 6 months. Visible FHB contamination on spikes was not observed in the field during any of the years concerned in this study.

Morphological characterization of fungal colonies

To determine the incidence of the different fungal genera, one hundred kernels were randomly chosen from each sub-sample, surface disinfected with 1% NaClO solution for 2 min, rinsed 4 times with running water, and dried for 12h under laminar flow hood. The use of running water cannot exclude a contamination introducing fungal spores and/or bacteria on the grains. Ten kernels were evenly placed into 10 separate Petri dishes (Ø 55 mm) containing potato dextrose agar, and incubated for 5 days at room temperature under conventional lighting without daynight cycle. Each colony was identified on cultural and morphological criteria to genus level among four genera, i.e. Alternaria, Epicoccum, Fusarium and Microdochium based on pigmentation of aerial mycelium and growth rate of mycelium, especially following the descriptions proposed by [START_REF] Pitt | Fungi and food spoilage, 2nd edn[END_REF]. The incidence was considered as the proportion of kernels infected by fungi of these four genera. During this long-term monitoring, we used the same morphological characterization of fungal colonies. However, we point out that the results of this study are relative to the method used. In particular, it is possible that some fungal species can be overgrown by colonies of other fungal species growing more rapidly, even if the incubation time was intended to limit this type of bias.

Incidence of different fungal genera on wheat grains in low-input and organic farming conditions The average incidence of each genus, or the absence of aerial mycelium development, is graphically represented per year and site, considering all bread wheat genotypes combined. This graphical representation constitutes a synthesis of the useful data for low-input and organic farming, i.e., 13 years × 3 sites × 5-16 genotypes per year for low-input farming conditions and 13 years × 2 sites × 5-8 genotypes per year in organic farming conditions. This graphical representation has not been subjected to statistical analysis due to the unbalanced design of the data, i.e., the lack of data in certain years on certain sites. The statistical analyses described hereafter were conducted on a suitable sub-dataset.

Sources of variability of the incidence of fungi on wheat grain per genus in low-input and organic farming conditions

For each fungal genus and farming condition separately, we assessed the effects of year, site, genotype and their interactions on the incidence of fungi on wheat grains using linear models.

First, some data were removed before applying statistical analyses, i.e., the genotypes tested in a single or two sites for low-input and a single site for organic farming conditions. The datasets focused on 7 years × 3 sites × 26 genotypes in low-input and 11 years × 2 sites × 13 genotypes in organic farming conditions (see Online Resource 2 for the list of genotypes).

Second, we implemented a full additive linear model including year, site, genotype and all second-order interactions, and we compared this full model with all other simpler additive models, including only a subset of the environmental variables. The comparison was based on the second-order Akaike information criterion (AICc; [START_REF] Sugiura | Further analysts of the data by akaike's information criterion and the finite corrections: Further analysts of the data by akaike's[END_REF], and considering all the explanatory variables as factors. We checked for homogeneity of variances (Bartlett test, visual inspection), normality of the residuals (Shapiro test, visual inspection) and influential points (Cook's distance). Sum to zero contrasts were applied prior to fitting models. We applied a transformation (log10+1, Box-Cox+1, Box and Cox 1964) or not on the response variable due to asymmetry in the data distribution depending on the fungal genus. Assumptions were met after transformation or not for the genera Fusarium, Epicoccum and Alternaria.

Third, the best model was retained, that corresponds to the most parsimonious model with the lowest AICc (AICbest). The best model was compared with the intercept-only model using AICc values, i.e., ΔAICc = AIC0 -AICbest in order to ensure the relevance of the best model (a value lower than 2 indicating that the best model was equivalent to the intercept-only model, [START_REF] Burnham | Multimodel inference: understanding AIC and BIC in model selection[END_REF]. The goodness of fit of the models was estimated using an adjusted R² (R²adj). We performed a type III ANOVA on the best model.

Concerning the genus Microdochium, no transformation was relevant to validate the assumptions. This result is due to the very strong asymmetry in the data distribution due to the absence of Microdochium observed in most cases (70% in low-input conditions, 52% in organic conditions), making the application of a linear model unsuitable. We have chosen not to present the results for the genus Microdochium in this study. First, this asymmetry requires the use of another type of model that would not be directly comparable with the results of the linear models (e.g. quasi-binomial generalized linear model). Then, there was non-independence of the four fungal genera in our study, the development of the aerial mycelium of a given fungal genus excluding the development of another genus (except in very rare exceptions), so showing the results of the four genera would lead to overinterpreting the data.

Statistical analyses were performed using R 4.0.1 (R Development Core Team 2022), its packages 'ggplot2' for graphical presentation, 'forecast' for Box-Cox transformation, 'AICcmodavg' for model selection based on AICc and 'car' for type III ANOVA.

Results

Incidence of different fungal genera on wheat grain in low-input and organic farming conditions

Data is based on 289 wheat samples collected under low-input and 155 wheat samples collected under organic farming conditions. Several main trends can be drawn from these long-term time monitoring. As stated above, it should be noted that the results are relative to the method used to characterize the fungal colonies.

First, the incidence of the genus Fusarium in bread wheat grains was generally low to nil, with a few rare exceptions (Fig. 1, Fig. 2). Indeed, the incidence of Fusarium was higher than 10% only in 2008 in organic conditions and in 2010 in low-input conditions, while lower than or equal to 1% in 61% of 'year × site' combinations in organic conditions and 64% of 'year × site' combinations in low-input conditions.

Second, the absence of aerial mycelium on bread wheat grain was not the most frequently observed case. Indeed, we observed a single combination with more than 50% of grains without fungal development in low-input conditions and no combination in organic farming conditions (Fig. 1, Fig. 2). On average, 89% of wheat grains were characterized by a fungal development in both low-input and organic conditions, with significant variability between 'year × site' combinations.

Third, in the vast majority of combinations considered, the genus Alternaria was dominant on wheat grains, i.e., in 89% of combinations in low-input conditions and 92% of combinations in organic conditions. The other genera showed a high incidence only in certain years and in certain sites. For instance, the genus Microdochium showed substantial development in Ile-de-France in 2016, while the genus Epicoccum was characterized by a higher incidence in 2012 (Fig. 1, Fig. 2). Considering the relevant sub-datasets and after confirmation that the prerequisites were met, a best model was identified for each case studied. This best model was always much more relevant than the intercept-only model, confirming that at least some explanatory variables are significantly explanatory for each case (Table 1). Moreover, R²adj was between 0.62 and 0.92, highlighting medium to high accuracies, depending on the model (Table 1). The statistical analyses carried out underlined that, for all the models tested, regardless of the fungal genus and the farming conditions, the year and the 'year × site' interaction had a large effect on the incidence of fungi on wheat grain. In contrast, the site effect was often non-significant or weakly significant, especially under organic conditions, and a low genotype effect was observed only for Alternaria in organic conditions (Table 1). It should be noted that a bias of the statistical analysis was related to the replacement of genotypes studied over time. Indeed, the long-term monitoring did not always characterize the same genotypes, this effect being then partially confused with the year effect, potentially explaining why the year effect is usually the only one significant in the analysis.

Table 1 Results of the best additive linear models indicating the effect of year, site, genotype and year × site second-order interaction on the incidence of the genera Fusarium, Epicoccum and Alternaria in both low-input and organic farming conditions considered separately. Box-Cox+1 and log10+1 refer to the Box-Cox transformation (transformation parameter λ = -1) and the base 10 logarithm after adding a constant (+1) to the raw values, respectively. AIC0 -AICbest corresponds to the difference of AICc between the best model and the interceptonly model. The goodness of fit of the models was estimated using an adjusted R² (R²adj). P-value estimation is based on F-test, *P < 0.05, ***P < 0.001, n.s., not significant. / indicates that the explanatory variable was not included in the best model. Degrees of freedom (df) vary between models applied on datasets in low-input and organic farming conditions

Discussion

This long-term study focusing on the incidence of four fungal genera on bread wheat grain over a 13-year period in several sites on hundred samples collected in both conventional low-input and organic farming conditions led to conclude two general trends.

First, the incidence of the genus Fusarium in bread wheat grains was generally low to nil, both in low-input and organic farming conditions, as observed in previous studies [START_REF] Lazzaro | Organic vs conventional farming: Differences in infection by mycotoxin-producing fungi on maize and wheat in Northern and Central Italy[END_REF][START_REF] Lenc | Fusarium head blight (FHB) and Fusarium populations in grain of winter wheat grown in different cultivation systems[END_REF]. In the vast majority of situations studied, the genus Alternaria was dominant on wheat grain, while the other genera studied were generally few or not observed, except in atypical years, sometimes only in one site. In our case, these particular results could be explained by meteorological conditions favourable to the development of fungal genera associated with maize as preceding crop (van der Velde et al. 2020), as also observed in other studies during part of the period considered here [START_REF] Giraud | Fusarium head blight and associated mycotoxin occurrence on winter wheat in Luxembourg in 2007/2008[END_REF][START_REF] Chandelier | Fusarium species and DON contamination associated with head blight in winter wheat over a 7-year period (2003-2009) in Belgium[END_REF]). It should be noted that in 2008 in low-input conditions in Brittany, i.e., 1km away from the organic trial characterized by a incidence near 50% of the genus Fusarium, considering a similar sowing period but faba bean as preceding crop, the incidence of Fusarium was under 5%.

Second, the main explanatory factors of the incidence of fungi on wheat grain, independently of the genus and the farming conditions, were the year and the 'year × site' interaction, while the site and the genotype were less explanatory. These results appear consistent with the trends observed on the graphical representations, suggesting the dominant influence of the environmental conditions of the year and the site on the relative development of the fungal genera on wheat grains. If certain conditions are known to favour the genus Fusarium, their combination have been rarely observed during our study period, except in 2008 in Brittany.

Indeed, in organic farming but also in low-input conditions, both the low level of nitrogen supply and a crop sequence including fewer cereals with maize rarely as the preceding crop, reduce the risk of FHBdevelopment (Dill-Macky and Jones 2000; [START_REF] Lemmens | The effect of nitrogen fertilization on Fusarium head blight development and deoxynivalenol contamination in wheat[END_REF][START_REF] Heier | Influence of N-fertilization and fungicide strategies on Fusarium head blight severity and mycotoxin content in winter wheat[END_REF][START_REF] Osborne | Epidemiology of Fusarium head blight on small-grain cereals[END_REF][START_REF] Bernhoft | Influence of agronomic and climatic factors on Fusarium infestation and mycotoxin contamination of cereals in Norway[END_REF]. Moreover, the main conducive weather conditions to FHB development are known, i.e., usually warm and moist conditions near and during the flowering period of winter wheat [START_REF] Osborne | Epidemiology of Fusarium head blight on small-grain cereals[END_REF]. However, the best predictors and window frames identified are partially different among studies [START_REF] Kriss | Relationship between yearly fluctuations in Fusarium head blight intensity and environmental variables: a window-pane analysis[END_REF][START_REF] Bernhoft | Influence of agronomic and climatic factors on Fusarium infestation and mycotoxin contamination of cereals in Norway[END_REF][START_REF] Landschoot | An empirical analysis of explanatory variables affecting Fusarium head blight infection and deoxynivalenol content in wheat[END_REF][START_REF] Shah | Predicting Fusarium head blight epidemics with weather-driven pre-and post-anthesis logistic regression models[END_REF]). In the absence of these conditions, the genus Alternaria predominated very frequently in our conditions.

Conclusions

Even though it should be noted that the results of this study are relative to the particular method used to characterize the fungal colonies, these results confirm the low to nil incidence of the genus Fusarium in most years in organic farming, just as in low-input farming conditions. These results are in accordance with previously published studies highlighting the low contamination level of mycotoxins, comparable or lower in organic than in conventional farming [START_REF] Brodal | Mycotoxins in organically versus conventionally produced cereal grains and some other crops in temperate regions[END_REF]. Moreover, from a breeder's point of view, this low incidence of FHB largely dependent on environmental conditions and agronomic choices [START_REF] Osborne | Epidemiology of Fusarium head blight on small-grain cereals[END_REF], makes phenotypic selection of more resistant genotypes in fields both more complex and less effective than for other diseases in northwestern Europe. Indeed, resistance testing appears to be more efficient under uniform and intermediate to high disease pressure, as usually achieved using artificial inoculation, and marker-assisted selection is widely put forward for FHB resistance [START_REF] Buerstmayr | Breeding for Fusarium head blight resistance in wheat -Progress and challenges[END_REF]). However, some additional morphological characteristics may also be considered in a breeding scheme, with taller plants characterized by anther extrusion being frequently more resistant to FHB [START_REF] Buerstmayr | Breeding for Fusarium head blight resistance in wheat -Progress and challenges[END_REF]. The results of the survey support that organic and low-input production of bread wheat does not represent a significant risk for the development of FHB in France. To go further in the characterization of Fusarium graminearum for future studies aimed at improving resistance to FHB, we strongly recommend the use of other methods [START_REF] Horevaj | A real-time qPCR assay to quantify Fusarium graminearum biomass in wheat kernels[END_REF][START_REF] Kumar | Real-time quantitative PCR based method for the quantification of fungal biomass to discriminate quantitative resistance in barley and wheat genotypes to Fusarium head blight[END_REF].

Fig. 1

 1 Fig. 1 Average incidence of the four different fungal genera on bread wheat grain in low-input farming between 2007 and 2019 at three sites are represented by the following shades of gray from bottom to top: Fusarium in black, Microdochium in dark gray, Epicoccum in gray, Alternaria in light gray, and absence of aerial mycelium observed in white. All bread wheat genotypes per year and site have been used to compute the average incidence, i.e., between 5 and 16 genotypes depending on the year considered. NA indicates an absence of data monitored concerning the particular 'year × site' combination. Bri = Brittany, IdF = Ile-de-France, P-C = Poitou-Charentes

Acknowledgements The authors acknowledge the teams from Institut national de la recherche agronomique (INRA, currently INRAE) experimental units of Auzeville, Lusignan, Le Moulon, Le Rheu, as well as the team of Agri-Obtentions for their involvement in these experiments.

Our work has been supported for 14 years by French Technical Institute for Organic farming (Laurence Fontaine ITAB) and InterBio Bretagne as part of the State-Region (Brittany) Project

Data availability

The datasets and metadata will be accessible on the publicly available Data INRAE repository at https://entrepot.recherche.data.gouv.fr/dataverse/inrae

Supplementary Information

Online Resource 1 Coordinates of the sampling locations in conventional low-input and organic farming systems Online Resource 2 List of bread wheat genotypes retained for the statistical analyses in conventional low-input and organic farming systems