A neutron-leakage spectrum model for on-the-fly rehomogenization of nodal cross sections

Abstract Modeling spectral effects due to core heterogeneity is one of the major challenges for current nodal analysis tools, whose accuracy is often deteriorated by cross-section homogenization errors. AREVA NP recently developed a spectral rehomogenization method that estimates the variation of the assembly-averaged neutron flux spectrum between environmental and infinite-lattice conditions using a modal synthesis. The effectiveness of this approach is tied to the evaluation of the spectrum of the neutron leakage from or into the assembly in the environment. In this paper, we propose a method for the leakage spectral distribution building upon Fick’s diffusion law. The neutron-exchange spectrum at a nodal interface is computed as a function of the gradient of the environmental flux spectrum, which is determined by the rehomogenization algorithm. This diffusive approach is applied to PWR benchmark problems exhibiting strong interassembly heterogeneity. We show that the method accurately reproduces the energy dependence of streaming effects, and that significant improvements in the input nodal cross sections, fission power and multiplication factor estimates are achieved at a low computational cost. The proposed model is compared with an alternative approach, that uses the fundamental-mode leakage spectrum obtained from the solution of the B1 equations. This second strategy is generally less accurate, and can only provide an adequate approximation of the environmental leakage in weakly heterogeneous systems.


Introduction
Routine calculations for reactor core design, monitoring and safety analyses are commonly performed with advanced nodaldiffusion methods on coarse meshes (Lawrence, 1986;Stacey, 2007). Fuel-assembly homogenization for the generation of fewgroup constants (nodal cross sections and discontinuity factors) is performed via heterogeneous transport calculations under the assumption of reflective boundary conditions at the assembly outer edges (Smith, 1986). However, this approximation can lose its validity when the assembly is simulated within the real environment (i.e., the reactor core). Here, streaming effects induced by internodal heterogeneity can cause significant deviations of the actual neutron flux distribution from the infinite-medium one used for spatial homogenization and energy collapsing of cross sections. Common examples in which the homogenization error can be highly penalizing are configurations with strong burnable absorbers and control rods; mixed oxide (MOX) assemblies sur-rounded by uranium oxide (UOX) assemblies; fresh-fuel assemblies facing depleted regions; and fuel bundles bordering reflector nodes. With these diverse layouts, the equivalence between the homogeneous nodal representation and the heterogeneous fine-mesh transport solution is only ensured if environmental (spatial and spectral) effects are modeled.
Several methods can be found in the reactor physics literature to correct single-assembly cross sections for spectral effects. Among them, we mention: empirical correlations taking into account local spectral interactions (Palmtag, 1997;Ban and Joo, 2016;Smith, 2017); the parameterization of nodal cross sections and discontinuity factors versus the current-to-flux ratio (and/or other albedo parameters) at the node outer surfaces (Rahnema and Nichita, 1997;Kim et al., 2017); high-order cross-section homogenization (Rahnema and McKinley, 2002); a spatial superposition technique of typical four-assembly configurations (Clarno and Adams, 2005); a recondensation method based on the Discrete Generalized Multigroup (DGM) energy expansion theory (Zhu and Forget, 2011); and a semi-heterogeneous transportembedded approach (Groenewald et al., 2017). The present work builds upon the spectral rehomogenization method developed at AREVA NP (Dall'Osso et al., 2010;Gamarino et al., 2017Gamarino et al., , 2018. In this approach, the variation of the neutron flux spectrum in the homogenized assembly between the environmental and infinite-medium conditions is estimated during the core nodal calculation via modal synthesis. The energy-condensation defects are computed on-the-fly and added to the few-group nodal cross sections interpolated from the standard parameterized tables. The performance of the method depends on two paramount points: (i) the set of basis and weighting functions employed for the modal expansion of the spectrum perturbation, and (ii) the definition of an accurate spectral distribution of the neutron leakage in the real environment. The former topic has been extensively treated in previous work (Gamarino et al., 2017(Gamarino et al., , 2018. In Gamarino et al. (2018), two modal approaches have been investigated. The first strategy uses analytical basis functions (Chebyshev polynomials of the first kind) and a physical mode in the fast group (i.e., the neutron fission-emission spectrum). The second approach is based on the Proper Orthogonal Decomposition (POD). It computes the optimal (in a least-squares sense) orthonormal basis functions for the space spanned by a set of snapshots of the reference spectrum perturbation. The two methods have been compared in terms of accuracy and computational efficiency. Several aspects of the rehomogenization method have been discussed, such as the implementation features, the impact of the approximations in the derivation of the algorithm, and the complementarity with other kinds of cross-section corrections (i.e., spatial rehomogenization and the critical-buckling spectrum correction).
In this paper, the methodology for the leakage spectral distribution is described. The leakage rate in a fuel assembly is dominated by two factors (Hebert, 2009): scattering anisotropy and interassembly neutron exchange. The former has an important effect in Pressurized Water Reactors (PWRs) due to the presence of hydrogen in the moderator, and is usually taken into account via transport corrections (such as the consistent B 1 and P 1 approximations) performed at the lattice-calculation level. The latter is inherently dependent on the core environment. The inaccurate results achieved with a flat-leakage approximation (i.e., considering the leakage spectral distribution uniform and equal to the coarse-group nodal estimate) highlighted the importance of finding a realistic energy shape for streaming effects (Gamarino et al., 2018). Hence, the aim of this work is to develop a model for the leakage spectrum. Two approaches are proposed and investigated. The first one is based on the application of Fick's diffusion law to the node-averaged environmental spectra estimated by the rehomogenization algorithm. We refer to it as diffusive-leakage model. The second one uses the homogenized-assembly criticalleakage spectrum from the fundamental-mode (B 1 ) calculation. The two strategies are tested on PWR assembly layouts characterized by significant heterogeneity. Both isothermal fresh-fuel conditions and configurations with depletion feedbacks are considered. Focus is given to the more promising diffusive-leakage approach.
This paper is structured as follows. In Section 2 the diffusive and fundamental-mode leakage methods are described. Section 3 shows numerical results for several PWR benchmark problems. In Section 4 we address various features of interest of the diffusive model. Concluding remarks and suggestions for future work follow in Section 5.

Description of the method
In this Section, the spectral rehomogenization method is briefly reviewed for the sake of completeness. The description of the two models for the leakage spectrum follows.

Review of spectral rehomogenization
The details about the derivation and the implementation features of the method can be found in Gamarino et al. (2018a).
For a generic homogenized node, the neutron continuousenergy balance equation in the environmental conditions can be written, within the coarse group G, as R t;G ðuÞU env;G ðuÞ þ L env;G ðuÞ ¼ The lethargy-like quantity u, bounded between 0 and 1, is defined as where E þ G and E À G denote the G-th group upper and lower energy boundaries. In Eq. (1), U env;G ðuÞ and L env;G ðuÞ represent the neutron spectrum and the leakage energy distribution, respectively. The remaining symbols have the conventional meaning (Stacey, 2007).
The assumption is made that the cross-section distributions depend only weakly on the environment (namely, R x;G ðuÞ % R 1 x;G ðuÞ for reaction type x). From now on, when referring to spectral functions we omit the argument u for the sake of lightness of the notation (i.e., In each of the N G coarse groups, the environmental spectrum is formulated as the sum of the reference distribution in the infinitemedium conditions ðu 1;G Þ and of the spectrum variation in the real environment ðdU G Þ: In Eq. (3), U G denotes the few-group node-averaged flux. The singleassembly spectrum u 1;G is normalized to unity, and dU G has zero average within G. The spectrum perturbation is expanded in terms of N Q G zero-averaged modal components Q G;i : Eqs.
The coefficients a G;i are solved for with a weighted-residual technique: after substitution of Eqs. (3) and (4), Eq. (1) is projected over a set of weighting functions W G;j (with j ¼ 1; . . . ; N Q G ) and integrated over u. The following N G N Q G Â N G N Q G linear system is obtained: In Eq. (5), the leakage projection coefficient c G;j is defined as whereas the remaining variables read as v G;j ¼ The reference (h R;x;G;j ) and variational (h V;x;G;i;j ) rehomogenization coefficients detailed in Eq. (7) depend on infinite-medium distributions and on the basis and weighting functions chosen for the modal synthesis of dU G . They are computed for each fuel-assembly type during the lattice calculation and stored in the cross-section libraries as additional homogenization parameters. During the online core calculation, the rehomogenization problem of Eq. (5) is solved for each node following a non-linear flux iteration. The U G and k eff estimates are taken as input from the previous, partially converged power iteration, and the rehomogenization coefficients are interpolated from the parameterized tables as a function of the local conditions. In conclusion, the spectral cross-section correction for reaction type x in a generic node is computed as follows: where the subscript 0 in h V;x;G;i;0 refers to the fact that W G;0 ðuÞ ¼ 1.

The leakage spectrum model
The leakage spectrum L env;G [Eqs.
(1) and (6)] is expressed as follows: where L G is the few-group node-averaged leakage and f L;G is a form function describing the leakage energy dependence. The distribution f L;G is normalized to unity so as to satisfy the condition Using Eq. (9), Eq. (6) becomes In the following, we formulate f L;G and h L;G;j for the two leakage models.

The diffusive-leakage method
We consider two adjacent nodes k and l separated by a surface DS along the generic direction x (Fig. 1). The two nodes have size Dx k and Dx l along x.
We apply the discrete (in space) Fick's diffusion law to compute the spectral distribution of the neutron current J S G through the surface DS: where D k G denotes the distribution in energy of the spatiallyaveraged diffusion coefficient in node k; U S env;G is the environmental spectrum at the interface between the two facing nodes, and U kc env;G is the environmental spectrum at the center of node k. A similar equation can be written for node l: We make the approximation that the spectrum at the center of a given node is equal to the node-averaged spectrum: The discrete formulation of Fick's law [Eqs. (12) and (13)] and Eq. (14) are based on a linear flux spatial distribution. This hypothesis is not consistent with the quartic polynomial expansion commonly adopted in advanced nodal codes. Because of the lack of information for a more rigorous spatial discretization of spectral distributions, we make the assumption that this approximation is acceptable within the range of accuracy of the proposed methodology.
As done for the cross-section distributions in Section 2.1, the dependence of the fine-energy diffusion coefficient on the environment is neglected: Continuity of the current distribution J S G is imposed by equating Eqs. (12) and (13). The following expression is found for the spectrum at the surface DS: where the quantityD m G ðuÞ reads as with the harmonic-averaged diffusion parameterD k;l G defined aŝ We refer toD k;l G as nodal-coupling diffusion coefficient. Moving to a more general multi-dimensional framework, the node-averaged leakage spectrum for the homogenized region k is determined applying Eq. (18) to all the interfaces with the surrounding nodes: The N nb normalization constants w k;m G in Eq. (21) are introduced to fulfill Eq. (10). They are computed imposing the preservation of the few-group leakage L k;m G through the interface between regions k and m: where L k;m G is defined in terms of the surface-averaged net current J k;m G at the interface: An estimate of J k;m G is known from the previous iteration of the nodal calculation. In Eq. (22), dividing by the node-averaged leakage L k G is required to scale f k L;G to unity. After introducing Eqs. (3) and (4), Eq. (21) reads as f k L;G ðuÞ ¼ The projection of Eq. (24) over the weighting functions W G;j (with j ¼ 1; . . . ; N Q G ) leads to the following definition of the j-th leakage rehomogenization coefficient for node k: In Eq. (25), a new type of rehomogenization parameter has been introduced for the nodal-coupling diffusion coefficient: The normalization condition of Eq. (22) results in the following expression for w k;m G : where, as in Eq. (8), we have used the fact that W G;0 is equal to unity to define the rehomogenization coefficients for j ¼ 0.
As observed in Eq. (26), the coefficients h R;D and h V ;D for a certain node k are not uniquely defined. This is because they also depend on the reference collapsing spectrum ðu m 1;G Þ and on the diffusioncoefficient distribution ðD m G Þ in the adjacent node m. The information on the former is carried by the coefficient h k;m Ã R;D;G;j , whereas the information on the latter is present in the coefficients h k Ã ;m R;D;G;j , h k;m Ã R;D;G;j and h k;m V;D;G;i;j . For a given fuel assembly, these ''mixed" rehomogenization parameters must be computed during the lattice calculation for each spectral interface (namely, for each dissimilar adjacent assembly). Nevertheless, the neighbor-bundle information is not easily achievable in the single-assembly simulation, because the lattice code has no knowledge of the bordering regions that the fuel assembly will see during its operating life in the reactor core. Although the cross-section generation procedure could be modified to add such a feature, this would demand to redefine the architecture of the lattice code. Moreover, the simulated unit assembly and its neighbor elements can experience different operating conditions and burn-up. Thus, the coefficients of Eq. (26) should be computed for several combinations of values of the state parameters in the adjacent assemblies, with the fuel exposure being the most relevant quantity to be sampled. In the light of the complex assembly-shuffling strategies adopted in modern core designs, the growth of the cross-section libraries caused by the storage of the rehomogenization parameters for the nodalcoupling diffusion coefficient would be significant. In conclusion, the formulation of the diffusive-leakage model presented above is not suitable for a practical integration in lattice-physics codes.
In order to overcome this issue, a variant of the method is proposed. It is based on the assumption that the diffusion-coefficient spectral distribution does not change significantly in the two adjacent assemblies: Using Eq. (28) again, the following expression holds for J S G : We emphasize that, although the approximate definition of Eq. (31) is not physical, it is justified by the fact that D G is almost not varying. The values of J S G computed for nodes k and l with Eq. (31) are equal and opposite. Hence, the continuity of the neutron-current spectrum at the interface DS is satisfied. Based on Eq. (31), the following approximation for f k L;G ensues: In Eq. (32), the spatial coefficient s k;m is given by After some algebraic manipulation, the leakage projection parameter becomes An analogous equation can be written for the same coefficients in the generic neighbor node m. Finally, the normalization constant The variables detailed in Eq. (35) are the standard rehomogenization parameters for the diffusion coefficient. They only depend on the infinite-medium neutron spectrum and diffusion-coefficient distribution in the given assembly. No information on fine-group quantities in the neighbor nodes is required. Therefore, they can be easily computed during the lattice calculation, in a similar manner as the rehomogenization coefficients for the cross sections and the fission spectrum (Eq. (7)).
Despite the heuristic connotation of Fick's law, the diffusive approach has a physical justification. This can be illustrated with an example. We consider a 3.1%-enriched UO 2 assembly with burnable-absorber (Pyrex) rods adjacent to a 1.8%-enriched UO 2 assembly. Both fuel bundles are at zero burn-up. For the assembly with poison, Fig. 2 shows the leakage form functions computed by Eqs. (21) and (32) (that is, the original and approximate formulations of the method) using the reference environmental flux spectra (i.e., U ref env;G ). The comparison with the reference environmental leakage reveals that the diffusive definition provides a very accurate estimate, and that the differences between the two formulations are negligible. In Sections 3 and 4 we discuss further the validity of the approximation of Eq. (28) and its effect on the calculation of the leakage spectrum.
From a numerical point of view, the diffusive approach translates into the dependence of the leakage parameter h L;G;j [Eqs.
(25) and (34)] on the modal-expansion coefficients a G;i , which are the unknowns of the rehomogenization algorithm. Therefore, a non-linearity is introduced. In addition, the spectral-correction problem is no longer local, because the spectrum-variation solution in a given node depends on the spectrum perturbation (i.e., on the coefficients a G;i ) in the neighbor nodes. More details about numerical aspects of the method are given in Section 4.1.

The fundamental-leakage approach
The second approach consists of using the fundamental-leakage spectrum determined at the single-assembly calculation level.
We make the approximation where f 1;k L;G is the leakage distribution making the infinite lattice critical (i.e., k 1 ¼ 1). This function is computed in most latticephysics codes, in which the unit assembly is simulated under critical conditions. In the absence of information on the exact operating conditions and on the materials surrounding the assembly, this assumption provides the most realistic representation of the critical core environment. Commonly, the critical-leakage calculation is based on the homogeneous fundamental-mode B 1 approximation. An exhaustive description of the corresponding theory can be found in Hebert (2009). In our work, we adopt the following formulation of f 1 L;G (the superscript k is omitted): where B 2 is the critical buckling (i.e., the buckling coefficient enforcing a multiplication factor equal to unity), and D G is the G-th group leakage-coefficient spectrum (function of B 2 ). Both quantities come from the solution of the homogeneous B 1 equations (Hebert, 2009).
Their product (D G B 2 ) is the critical-leakage cross-section distribution. In Eq. (38), u B 2 1;G denotes the B 2 -corrected infinite-medium spectrum, which has the same shape in energy as the fundamental mode computed by the B 1 model. The normalization of f 1 L;G to unity satisfies Eq. (10). After substitution of Eqs. (37) and (38) into Eq. (11), the leakage projection coefficient for a generic node reads as With this approach, h L;G;j can be computed on the basis only of lattice information. Therefore, its calculation is performed directly during the single-assembly simulation, as for the other rehomogenization parameters (Eq. (7)). No complexity is added to the online solution of the spectral rehomogenization problem. Despite its simplicity, this method presents some significant limitations. Even if the B 1 model provides the best possible representation of the critical lattice surrounding the assembly, the infinite-medium shape formulated in Eq. (38) might not capture the streaming effects occurring in the real environment in the presence of strong interassembly heterogeneity. Moreover, the consistency of the B 1 spectrum correction fades when non-critical conditions are simulated, such as reactor core transients and subcritical states during reactor start-up or power outage. In these situations (Dall'Osso, 2015a,b;Demaziere, 2016), the B 2 -corrected spectrum and the fundamental-leakage distribution can differ from those in the non-critical core environment even in homogeneous systems (i.e., in the absence of streaming effects). Another drawback of this approach is its lack of generality, because it can only be applied if the cross-section libraries are built with the fundamental-buckling paradigm.

Numerical results
In this Section, the methodology is applied to two-group nodal simulations of several PWR examples. The analysis is made on colorset configurations, consisting of four-assembly sets with reflective boundary conditions at the assembly center-lines. In the first part, the diffusive-leakage model is validated. In the second part, the fundamental-leakage approach is tested and the two strategies are compared.

Validation of the diffusive-leakage model
Reactor configurations at initial-core isothermal conditions (i.e., without thermal-hydraulic feedbacks and fuel depletion) are first addressed. We consider the following benchmark problems: a UO 2 colorset with burnable-poison (Pyrex) rods (Example 1); a UO 2 colorset hosting silver-indium-cadmium (AIC) control rods (Example 2); a UO 2 /MOX colorset (Example 3); a UO 2 colorset with gadolinium fuel pins (Example 4). In Gamarino et al. (2018), the authors have used Examples 1, 2 and 3 to validate the modal synthesis of the spectrum variation. For the above test cases, nodal simulations are run with BRISINGR, a Delft University of Technology in-house-developed code. The solution strategy in BRISINGR is based on a conventional non-linear coupling between a Coarse Mesh Finite Difference (CMFD) solver and a Nodal Expansion Method (NEM) with fourth-order polynomial synthesis of the two-group intra-nodal transverse-integrated flux. Two-group homogenization parameters are computed with the continuousenergy Monte Carlo neutron transport code SERPENT (Leppanen et al., 2015). The details about their calculation can be found in Gamarino et al. (2018). The diffusion coefficients are formulated with the Cumulative Migration Method (CMM) (Liu et al., 2016). For the sake of generality, in this part of the analysis the singleassembly cross sections are generated without the criticalbuckling (B 2 ) correction. This approach is of particular interest in the light of our previous findings (Gamarino et al., 2018), showing that rehomogenization can also reproduce spectral effects due to different multiplicative properties in the core environment and in the infinite-medium lattice.
As further validation of the methodology, we also analyze a test case with fuel depletion (Example 5). This benchmark problem consists of a UO 2 colorset with Pyrex rods and is modeled with ARTEMIS (Hobson et al., 2013), the core simulator of AREVA NP's code platform ARCADIA (Curca-Tivig et al., 2007). The crosssection libraries used by ARTEMIS are generated with the deterministic lattice transport code APOLLO2-A (Martinolli et al., 2010).
In all the example problems, a nodalization of 2 Â 2 nodes per assembly is chosen. The values of the main state parameters correspond to standard hot full power (T fuel ¼ 846 K, T mod ¼ 582 K, p = 158 bar). For each benchmark problem, the results of the following calculations are presented: with infinite-medium cross sections (a); with cross sections corrected by the reference spectral defect (b); with spectral rehomogenization of cross sections using the ref- Rehomogenization is applied with Galerkin projection of Eq. (1) and the two kinds of basis functions investigated in Gamarino et al. (2018): Chebyshev polynomials of the first kind, in combination with a physically justified mode (the neutron emission spectrum from fission) in the fast group; proper orthonormal modes, computed via Singular Value Decomposition (SVD) of a set of snapshots of the reference spectrum variation in Examples 1, 2 and 3.
The snapshots for the calculation of the POD modes have been obtained parameterizing the fuel-assembly composition, namely, sampling several values of the 235 U enrichment, burnable-poison concentration, plutonium content in the MOX bundle, number and type (AIC, B 4 C) of control elements in the rodded configuration. Unless stated otherwise, in the examples which follow we use the polynomial synthesis in the fast group and the POD basis in the thermal group (hybrid approach). This choice is made in view of an application of the methodology to industrial reactor calculations. As found in Gamarino et al. (2017Gamarino et al. ( , 2018, the thermal spectrum variation exhibits a weak dependence on the type of spectral interface and on the local conditions (such as the fuel exposure). Therefore, the proper orthonormal modes computed by a limited set of snapshots can successfully synthesize it in several core configurations, even if samples of their solution have not been included in the snapshot array. On the other hand, in the fast group the spectrum deformation strongly depends on the assemblyinterface type and on the burn-up. Hence, the accuracy of the POD-based rehomogenization is tied to an extensive sampling of heterogeneous assembly configurations and fuel evolutions. Since at the current stage of development of the methodology a universal POD basis has not been achieved yet for the fast group, the analytical approach is deemed to have a more general validity in this energy range.
The rehomogenization coefficients [Eqs. (7), (25), (34) and (39)] are computed with a 281-group energy structure. The numbers of fine groups g used in the fast and thermal coarse groups are 247 and 34, respectively. The upper boundary of the fast group is (2)]. The lower boundary of the thermal group is E À 2 ¼ 1:1 Á 10 À10 MeV. The thermal cut-off energy is where dU ref G denotes the reference spectrum variation from the lattice code. With this choice, the spatial effects of the environment are not taken into account (Dall'Osso, 2014;Gamarino et al., 2016). In calculation c, the environmental leakage is computed with Eq. (38), using the assembly leakage cross-section distribution and flux spectrum obtained from a 281-group transport simulation of the whole colorset. Calculation c provides the reference solution to assess the performance of the leakage model, whereas calculation b provides the reference solution for the rehomogenization method as a whole.

Example 1 -UO 2 colorset with Pyrex rods
The colorset is made of four 17 Â 17 PWR fuel assemblies of fresh UO 2 with two different compositions: 1.8% enrichment, 3.1% enrichment and 16 rods containing burnable poison. The absorber elements are made of borosilicate glass (Pyrex). The colorset and assembly layouts are depicted in Fig. 3. The concentration of diluted boron in the moderator is 1465 ppm and corresponds to the critical value (i.e., k ref eff ¼ 1:0). The reference (normalized) total fission power is 0.92 in the 1.8%-enriched assembly and 1.08 in the 3.1%-enriched assembly with Pyrex. Fig. 4 shows the leakage spectrum, computed by rehomogenization with the diffusive model, in the assembly without burnable absorber. The curves are normalized to the few-group assembly-averaged leakage from the nodal calculation. Units are in neutrons/cm 3 /s. The reference environmental leakage from the transport calculation is accurately reproduced in the fast group. Minor deviations occur only in the high-energy peak range (for u 2 ½0:7; 0:95, i.e., approximately between 0.12 and 8.2 MeV) and in the epithermal region (for u < 0:1, i.e., E < 3:5 eV). In the thermal group the result is also satisfactory, even if a slight shift of the bell-shaped curve towards higher values of u is observed. The variations between the original and approximate definitions of the leakage function are negligible. Fig. 5 depicts the spectrum variation in the two assemblies. The percent relative perturbation is calculated with respect to the assembly-averaged two-group nodal flux. The curves computed with the diffusive approach have accuracy comparable to those ensuing from the reference-leakage input. The slight overestimation (in absolute value) of the reference perturbation in the epithermal region is a consequence of the aforementioned discrepancy in the leakage spectrum predicted by the method. The shift found in the computed thermal leakage has no appreciable effect on the spectrum deformation.
The errors in the nodal cross sections are reported in Tables 1 and 2 for the two assemblies. The corrections computed with the diffusive-leakage model reproduce very accurately those obtained with the reference leakage. A slight miscorrection is only found in fast-to-thermal scattering with the original formulation (calc. d). Table 3 shows the errors in the multiplication factor ðDk eff Þ, fewgroup nodal flux ðD U G Þ and nodal fission power ðDP fiss Þ. The two values reported for D U G refer to the fast and thermal groups, respectively. The value of DP fiss out of parentheses refers to the total power, and the two values within parentheses correspond to the fast-and thermal-group power, respectively. Also for these parameters, the deviations of simulations d and e are very close to those of the reference calculations (b and c). The error in fission power drops to zero. The residual errors in k eff and in the cross sections are due to the spatial effects of the environment, that are not taken into account by the method. Table 3 also indicates the number of non-linear power iterations ðN iter Þ for the convergence of the eigenvalue calculation. We used a tolerance iter ¼ 10 À5 for the relative variation of the k eff         3.1.2. Example 2 -UO 2 colorset with AIC control rods The colorset is composed of four 17 Â 17 UO 2 assemblies with 1.8% enrichment (Fig. 6). A bank of 24 AIC (silver-indiumcadmium) control rods is inserted in two of them. The moderator has no soluble boron ðC B 10 ¼ 0 ppmÞ. The reference multiplication factor is 0.98860, and the reference fission power is 1.22 in the unrodded assembly and 0.78 in the rodded one.
Figs. 7 and 8 show the leakage distribution and the spectrum variation determined by rehomogenization. In the thermal group, for both quantities the computed curves almost overlap with the reference-leakage ones. In the fast group, the leakage prediction is very precise in the epithermal range, whereas some inaccuracy arises in the high-energy region (for u 2 ½0:7; 1:0, that is, between 0.12 and 19.6 MeV). This causes a shift of the calculated spectrumperturbation peak towards higher values of u, as observed in Fig. 8. The fast-group spectrum deformation computed with the POD modes is plotted in Fig. 9 for the unrodded assembly (only the result of the approximate version of the leakage model is shown). Compared to the polynomial approach, the outcome is more accurate in the epithermal region, with the resonance spikes being fitted precisely.
The errors in the nodal cross sections, multiplication factor, nodal flux and fission power are shown in Tables 4-6. Since the results achieved with the two formulations of the diffusiveleakage model are equivalent, they are only reported for the approximate one. With the diffusive approach, the deviations in the thermal cross sections are very close to those found with calculation b. In the fast group, for all the reaction types the computed corrections go in the right direction and approach the reference ones. The corrections on k eff , the nodal flux and the fission power are also in good agreement with the reference values. The errors in the last two quantities are significantly lower than those found with infinite-medium homogenization parameters.
The error in the rodded-assembly R a;2 increases when the reference dR a;2 is added to the infinite-medium value (calc. b). This is due to the exclusion of spatial effects (Gamarino et al., 2018). Calculation c (i.e., with the reference leakage spectrum) somewhat  deviates from simulation b in the errors in fast-group absorption cross section and power. This is because the computed dU 1 does not capture the resonance spikes in the interval [0.1,0.2] (corresponding to E 2 ½3:5 eV;20 eV). Hence, an overcorrection on R a;1 occurs. In the simulation with the diffusive leakage and Chebyshev modes, the resonance peaks are better fitted. This partially compensates for the inaccuracies in the prediction of the global behavior, whose magnitude is slightly underestimated in the epithermal range. The small deviation in k eff in the calculation without rehomogenization is due to fortuitous error compensation.

Example 3 -UO 2 /MOX colorset
The third colorset, which is shown in Fig. 10, consists of two 18 Â 18 UO 2 and MOX assemblies. The UO 2 assemblies have 2.1%   enrichment. The MOX assemblies are made of three fuel-pin types differing in plutonium content and 235 U enrichment. The concentration of diluted boron in the moderator is 2907 ppm. The reference multiplication factor is 1.00194, and the reference fission power is 0.86 in the UO 2 assembly and 1.14 in the MOX bundle.
Figs. 11 and 12 depict the leakage distribution and the spectrum perturbation estimated by rehomogenization with the hybrid and POD modal approaches. Also in this case, the results are only shown for the approximate formulation of the leakage model. The computed fast-group distributions suffer from inaccuracy in the higher part of the energy domain (u > 0.85, i.e., E > 1.5 MeV). Here, the bulge-shaped outline featured by the spectrum variation is not reproduced by the polynomial approach, whereas it is amplified by the POD-based one. The result is instead satisfactory in the epithermal region. Tables 7-9 show the errors in the nodal cross sections and in the integral parameters. As in the previous examples, the performance of the method is excellent in the thermal group. In the fast group, all the cross-section corrections go in the right direction. The simulation with the POD modes (calc. e2) reproduces the reference dR a;1 almost exactly, whereas both calculations c and e1 overestimate the correction, especially in the MOX assembly. This difference depends on the reconstruction of the epithermal resonances, that is achieved to a high level of accuracy only with the POD basis. Due to the aforementioned flaws in the prediction of dU 1 at high energies, the three calculations overcorrect significantly the production cross section mR f ;1 and, as a result, the fast-group nodal power. Nevertheless, due to the small contribution of the latter and to the improvement in the thermal-power estimate, a considerable reduction of the error in the total power is found. The correction on k eff is overestimated (calc. e2) or mispredicted (calc. e1), but the errors remain low.

Example 4 -UO 2 colorset with gadolinium rods
In this example (Fig. 13), the checkerboard layout consists of two 17 Â 17 UO 2 assemblies with 1.8% enrichment and two 17 Â 17 UO 2 assemblies with 3.9% enrichment and 12 fuel rods containing gadolinium (Gd). The pins with burnable poison are located at the periphery of the assemblies and have 0.25% 235 U enrichment and 8% Gd enrichment. The concentration of boron in the moderator is 1830 ppm, and the reference multiplication factor is 1.00303. The reference fission power is 0.82 in the 1.8%-enriched assembly and 1.18 in the 3.9%-enriched one.
Figs. 14 and 15 show the leakage distribution and the spectrum variation. In the fast group, the spectrum change from the diffusive approach exhibits a tilt with respect to the curve obtained with the environmental-leakage input, and overestimates (in absolute value) the reference deformation in the epithermal region. The deviations in the computed dU 2 are due to spatial effects, and can be justified as follows. In the assembly with gadolinium the flux spatial variation is positive (and up to 15%) at the periphery, where the fuel pins with burnable poison are located and neutron absorption is higher. Hence, the global spatial correction dR spat a;2 in the node is positive. As the hardening effect of the spatial term is not accounted for by spectral rehomogenization, the method predicts a softer spectrum (that is, the amplitude of dU 2 is overestimated in the intermediate region of the thermal domain and in its upper part).  Table 9 Example 3: number of power iterations and errors in the multiplication factor, nodal flux and fission power. Tables 10-12 report the numerical errors. Simulations b and c differ most clearly in their prediction of dmR f ;2 in the assembly with Gd rods and of dR a;1 in both fuel bundles. In the thermal group, the corrections computed with the diffusive model match those of calculation c. The cross sections R a;1 and R s;1!2 are overcorrected due to the overestimation of dU 1 in the epithermal range. The correction on mR f ;1 is larger than the reference value in the assembly with poison, whereas it goes in the wrong direction in the low-enriched bundle. These inaccuracies are due to the misprediction of the spectrum change in the range [0.95,1.0] (i.e., E 2 [8.2 MeV,19.6 MeV]), where a non-zero dU 1 is computed. The errors in k eff and fission power do not decrease significantly if only spectral effects are accounted for.

Example 5 (depletion feedbacks)
The colorset is composed of three 1.8%-enriched UO 2 assemblies and a 3.1%-enriched UO 2 assembly hosting 16 Pyrex rods. The composition and the internal layout of the fuel bundles are the same as   those displayed in Fig. 3(b) and (c). The fuel elements are burnt at a power volumetric density of 107.03 kW/l until an average colorset exposure of 12.0 GWd/t (corresponding to about 303 days). The depletion is performed with 50 burn-up steps of gradually increasing size. The values of the state parameters are kept constant during the evolution. The diluted-boron concentration (1000 ppm) is chosen so as to make the configuration critical during the first part of the depletion (Fig. 16).
In this case, the two-group cross-section libraries are generated with the critical-buckling correction, which is the default option in the lattice code APOLLO2-A. Rehomogenization is applied with the Chebyshev basis functions also in the thermal group, and with the approximate variant of the diffusive-leakage method. In APOLLO2-A the rehomogenization coefficients are parameterized only versus burn-up. They are computed at predetermined nominal values of the fuel temperature, moderator temperature and density, diluted-boron and xenon concentrations. During the nodal calculation, they are updated to account for the differences between the actual values of the above state parameters in the node and the nominal ones. This choice is made to minimize the memory requirement for the storage of the additional homogenization entries. The algorithm developed for the update estimates the variation of the infinite-lattice condensation spectrum between a nominal and a perturbed state, using an approach similar to that described in Section 2.1. It requires to compute and store the isotopic rehomogenization coefficients of water, soluble boron and xenon in the nominal conditions. The full details of the update methodology will be given in future publications. In the framework of the present work, we verified that the error introduced using updated rehomogenization coefficients (instead of computing them at the exact local conditions) is negligible or small. Fig. 17 shows the spectrum variation in the 1.8%-enriched assembly next to the heterogeneous bundle (i.e., with Pyrex rods) at the beginning and at the end of the depletion. In the fast group, the reference curve exhibits a change of sign and a significant deformation with burn-up, especially at high energies. As observed for Example 3, rehomogenization succeeds in predicting the average behavior of the distribution in the epithermal and intermediate regions of the spectrum, but it fails to reproduce rigorously its strongly varying outline in the upper part of the energy domain. The comparison with the reference-leakage-input curve reveals that the leakage spectrum is accurately estimated by the diffusive model, and that the above inaccuracy is due to the inherent limitations of a polynomial synthesis in the fast group. In the thermal range, neither the shape nor the magnitude of the spectrum perturbation changes appreciably with the fuel exposure, and the reconstruction remains accurate throughout the evolution. Figs. 18 and 19 depict the errors in the absorption and production cross sections as a function of burn-up for the fuel assembly with Pyrex rods and the poison-free bundle next to it. In the plots the zero-error bar is highlighted. The corrections computed with the diffusive model are in good agreement with those obtained with the reference-leakage input. A significant overcorrection is only found for fast absorption in the assembly with burnable absorber. For both reaction types a considerable improvement is achieved in the thermal group of the heterogeneous assembly (above all in mR f ;2 ) and in the fast group of the assembly without Pyrex. In the poison-free bundle, errors in R a;2 are significantly reduced in the first part of the depletion. However, they increase with burn-up and ultimately overcome in magnitude the homogenization defect, which slowly decreases with the fuel exposure. No gain in accuracy is found for fast absorption in the heterogeneous assembly.
The errors in the nodal flux are shown in Fig. 20 for the above two assemblies and in Fig. 21 for the 1.8%-enriched assembly next to an assembly of the same type. The improvement produced by  rehomogenization is evident in the thermal group, where in the absence of spectral corrections the deviations increase significantly with burn-up (up to 2% in the assembly with Pyrex rods and about 1% in the remaining two bundles). With the diffusive-leakagebased rehomogenization, the errors are bounded below 0.5% in the dissimilar bordering assemblies and 0.2% in the third bundle. Furthermore, they do not change significantly with the fuel exposure, whereas in the calculation with infinite-medium cross sections they exhibit a monotonically increasing behavior during the second part of the depletion. Fig. 22 shows the evolution of the error in the nodal fission power for the unlike adjacent assemblies. The benefits of rehomog-enization are apparent, especially in the first part of the depletion. The behavior of the curves can be interpreted as follows. When infinite-medium cross sections are used, the power is undervalued in the more reactive assembly (i.e., the 3.1%-enriched one). This is due to the underestimation of its thermal production cross section (Fig. 19). Therefore, the fuel initially burns less and loses reactivity more slowly, which goes in the direction of an increase in the power with burn-up. The opposite occurs for the less reactive assembly type (1.8%-enriched), in which the power is overestimated. A consequence of this evolution is that the power deviations tend to decrease with burn-up in the three assemblies. As shown in Fig. 23, the error in the multiplication factor (which is  initially negative) similarly becomes lower for increasing values of the fuel exposure, and approaches a constant value. When rehomogenization is applied, the deviations in the power are considerably lower, especially in the assembly with Pyrex rods at the beginning of the depletion. However, with the error pattern introduced by the spectral corrections the aforementioned self-healing effect vanishes. This might be the cause (or one of the causes) of the monotonically increasing behavior of the deviations in the power and k eff (Fig. 23) when rehomogenization is applied. Another possible source of inaccuracy in the depletion is that the spectral corrections are only computed for the macroscopic cross sections. Few-group microscopic cross sections are not rehomogenized. Therefore, the solution of Bateman's equations for the depletion of fissile isotopes and of the burnable poison benefits from rehomogenization only in part (namely, via the improved accuracy in the few-group nodal flux, as observed in Figs. 20 and 21). This source of error could be removed introducing isotopic rehomogenization coefficients to correct the microscopic cross sections. For the generic nuclide c (and reaction type x), the reference and variational isotopic parameters can be defined as h R;x;c;G;j ¼ Rehomogenization of microscopic cross sections will be addressed in future work.

Comparison with the fundamental-leakage approach
In this part of the analysis, the fundamental-leakage approach is investigated and compared to the diffusive method for some test cases without feedbacks. The nodal calculations are performed with BRISINGR. The infinite-medium cross sections, discontinuity factors and rehomogenization coefficients are computed in APOLLO2-A. This choice has been made to avoid a computationally demanding fine-group B 1 spectrum calculation in SERPENT. The diffusion coefficients are obtained from the homogeneous B 1 model. This approach makes the assumption that the leakage coefficient D G defining the critical leakage (Eq. (38)) can be used as diffusion coefficient (Hebert, 2009). The results are briefly presented for Examples 1 and 2, which are now simulated using singleassembly input data generated with the critical-buckling procedure. An additional test case (a supercritical UO 3 /MOX colorset) is also considered.
For the assembly of Example 1 hosting Pyrex rods, Fig. 24 compares the fundamental-mode leakage computed in APOLLO2-A and the leakage predicted by rehomogenization with the diffusive model. The spectrum perturbation determined with the two approaches is also depicted. In the thermal group, the criticalleakage spectrum significantly overestimates the reference (in amplitude) for u 2 ½0:85; 1:0 ðE 2 ½0:15 eV;0:625 eVÞ, and underestimates it in the remaining part of the domain (u 2 ½0:5; 0:85, corresponding to E 2 ½6 meV;0:15 eV). As a consequence, the magnitude of the spectrum change is underestimated in the two lethargy ranges. In the epithermal region of the fast group (u 2 ½0; 0:3, E 2 ½0:625 eV;110 eV), the environmental leakage is negative (i.e., there is an incoming flow of neutrons), whereas the fundamental-mode leakage is positive. Therefore, the dU 1 computed with the latter deviates significantly from the reference in this lethargy range, and eventually has opposite sign  at the border with the thermal group. As in this region the finegroup cross sections are higher, an error in the prediction of the spectrum variation has more weight. In the high-energy region (u 2 ½0:75; 0:85, that is, E 2 ½0:27 MeV;1:5 MeV), the underestimation of the environmental leakage causes a considerable overprediction of dU 1 (in absolute value). Analogous results are found for the assembly without Pyrex rods. The effect of the mispredictions in the fast group can be observed in Table 13, showing the deviations in the integral parameters and in the main nodal cross sections (errors are in pcm for the multiplication factor and in percent for the fission power and cross sections). With the fundamental-leakage approach (calculation f), R a;1 is significantly overcorrected, whereas the corrections on mR f ;1 go in the wrong direction. For both reaction types, the errors become higher than in the calculation without rehomogenization. In the thermal group the corrections have the right sign, but their magnitude is underestimated. The error in the fission power increases notably. Also in this case, with the diffusive-leakage approach the cross-section corrections are accurately predicted and the errors in k eff and in the fission power match the reference ones.
The same analysis is performed for the colorset of Example 2. Fig. 25 depicts the leakage distribution and the spectrum variation in the rodded assembly. Compared to the case without criticalbuckling correction (Fig. 8), the change in the shape (and sign) of the fast-group spectrum deformation is apparent. The reconstruction of the perturbation with the fundamental-leakage approach still lacks accuracy. In the thermal group, similar conclusions can be drawn as for the previous test case. The deviations in the nodal cross sections and in the integral parameters are in Table 14.
Calculation e accurately corrects R a;1 , which is instead largely undercorrected with simulation f. The error in k eff is reduced by rehomogenization (especially with the diffusive-leakage model), whereas the impact on the fission power is small (most of the error is due to spatial effects).
In the examples considered so far, the errors in the fast-group infinite-medium cross sections are mostly due to spectral effects, rather than spatial ones. Hence, a last benchmark problem has been tailored to achieve very high errors in the fast-group cross sections. In this way, we can better evaluate the capability of rehomogenization to correct them. The example consists of a UO 3 /MOX colorset. Its layout is the same as that illustrated in Fig. 10(a), with one of the two MOX assemblies replaced by a UO 3 assembly. The three UO 3 assemblies have 3.5% 235 U enrichment, whereas the MOX assembly has 8.0% 239 Pu enrichment. The internal loading of the fuel bundles corresponds to that of Fig. 10(b) and (c). The diluted-boron concentration is 700 ppm, and the reference multiplication factor is 1.26257. The reference fission power is 0.87 in the MOX assembly, 1.03 in the UO 3 assembly next to it and 1.10 in the third UO 3 bundle. Figures 26 and 27 show the spectrum variation in the adjacent UO 3 and MOX assemblies (heterogeneous interface) and in the UO 3 assembly next to another UO 3 bundle (homogeneous interface). Compared to the previous test cases, the fundamentalleakage spectrum provides a better approximation in the fast group. In particular, the prediction is reasonably accurate for the UO 3 assembly next to the MOX assembly. However, in the latter the perturbation is considerably underestimated in the epithermal region. In the UO 3 assembly bordering an identical  assembly, the two leakage models produce a similarly accurate outcome. In this case, due to the low interassembly heterogeneity the fundamental-leakage approach provides a satisfactory approximation of the environmental leakage. In the thermal group, the infinite-medium leakage strategy is still inadequate, especially in the MOX assembly. Table 15 shows the deviations in the nodal cross sections and integral parameters. Their magnitude is relevant for the fast-group cross sections (up to 4% for absorption in the MOX assembly). A considerable part of the errors is due to the use of a B 2 -corrected collapsing spectrum in an environment strongly deviating from criticality. The benefits of rehomogenization in combination with the diffusiveleakage model are apparent. The errors in group-one cross sections are significantly reduced, and become very close to those of calculation c. An improvement in the power prediction is also found. With the fundamental-leakage approach the errors also decrease, but to a smaller extent. The deficiencies in the calculation of the thermal corrections cause a strong overprediction of the power in the MOX assembly. If no spectral correction is performed, the errors in the thermal flux (not shown in Table 15) are relevant: À6.67% in the MOX assembly, À4.88% in the adjacent UO 3 assembly and À4.92% in the remaining UO 3 bundle. These values are respectively reduced to À1.88%, 0.17% and 0.03% with calculation e, and to À2.88%, 0.87% and 0.15% with calculation f. In the second assembly, with calculation f the somewhat high residual error in U 2 deteriorates the estimate of the total power,   even if the deviations in R a;2 and mR f ;2 are close to the ones of calculation c due to fortuitous error cancellation. The variation of k eff in calculation c goes in the opposite direction with respect to calculation b. This is found, to a lesser extent, also for the fission power in the MOX assembly. Such outcome is again a side effect of the use of smoothly-varying basis functions, which cannot reproduce in any way the fine details of the spectrum perturbation (in particular, the resonance spikes). As final remark, in all the test cases the fundamental-leakage approach exhibits a somewhat faster convergence than the diffusive one (namely, the increase in the number of flux iterations is well below a factor of 2). This can be explained by the absence of non-linearity.

Discussion
In this Section we discuss some aspects of the diffusive-leakage method. Focus is given to its numerical features and computational efficiency and to the validity of the assumption underpinning the approximate formulation of the model (Eq. (28)). We also address the influence of the fine-group diffusion-coefficient definition on the accuracy of the leakage spectrum prediction.

Numerical aspects and implementation features
As mentioned in Section 2.2.1, the diffusive-leakage method is non-linear and non-local. Non-linearity is tackled computing the   The spectral correction problem cannot be solved independently for each node, as in the case of the fundamental-leakage approach.
This requires to store the coefficients a G;i for all the nodes of the system. However, with the Galerkin-based modal synthesis this must be done irrespective of the leakage model to dampen numerical oscillations observed in the convergence process (Gamarino et al., 2018). Therefore, non-local features of the diffusive approach add no further complexity to the rehomogenization problem.
When the normalization factor w k;m G is computed for the leakage spectrum at the interface between nodes k and m, numerical instability can arise if the denominator of Eq. (27) or Eq. (36) has a value close to zero. This can occur, for instance, in adjacent fresh-fuel nodes having the same composition. In order to avoid undesired convergence issues in the nodal simulation, the diffusive approach is only applied if the following approximate threshold condition is met: where tol is a given tolerance and U avg;km G is the arithmetic average of the two-group flux in nodes k and m. Otherwise, no action is performed (i.e., the assumption is made that there is no leakage through the given interface) or the fundamental-leakage shape is used for that surface, if available from the lattice calculation. Values of tol in the range 0.1%-1.0% proved to be a suitable choice. A similar threshold condition should be applied to the node-averaged leakage L k G , that divides J k;m G in Eqs. (27) and (36). However, this second control is not necessary if we use the notation of Eq. (20) with normalization to J k;m G , instead of the notation of Eqs. (9) and (21) with normalization to J k;m G = L k G . Due to the non-linearity of the diffusive approach, the spectrum perturbation and the leakage distribution computed by rehomogenization are intimately coupled. In order to reduce oscillations in the convergence of the calculation, numerical damping is performed on the leakage coefficients at each assembly interface.
We define the leakage parameter q k;m G;j for the edge of node k bordering with m (we take as example the approximate formulation): where q k;m ðnÀ1Þ 0 G;j denotes the estimate of Eq. (44) from the previous rehomogenization iteration, and h is the damping factor (taken as 0.5). Under-relaxation has been found to strongly increase the stability of the computation. Without it, for most test cases a significantly lower rate of convergence is observed. The drawback of this procedure is that the coefficients q k;m ðnÞ 0 G;j must be saved in memory for each internodal surface. For a three-dimensional two-group simulation with four basis functions per group, the memory requirement is 48 coefficients per node. The results of Section 3 showed that the computational cost of spectral rehomogenization combined with the diffusive method is reasonably low: for most of the test cases, the increase in the number of non-linear flux iterations is below a factor of 2.

On the approximate formulation of the method
As highlighted in Section 3, the differences between the original and approximate formulations of the diffusive method are negligible. This confirms that the approximation of Eq. (28) is generally acceptable. For example, Fig. 28 depicts the fine-group diffusion coefficient in the two assemblies of Example 1. The deviations between the two distributions are negligible in the fast group, whereas they are more substantial in the thermal group. However, the rehomogenization method is less sensitive to variations of the thermal leakage, which has smaller magnitude.
The behavior observed in Fig. 28 can be justified considering that the migration of fast neutrons is not significantly affected by differences in the enrichment or by the presence of thermal absorbers (such as burnable poison and control rods). These are only perceived when neutrons are slowed down to thermal energies. A similar outcome has been also found when comparing the diffusion-coefficient distributions of adjacent assemblies with significantly different fuel composition, such as the UO 2 and MOX assemblies of Example 3.
To conclude, we verified that environmental effects on the fineenergy diffusion coefficients are negligible, and that the performance of the method is not influenced by the use of the infinite-

Impact of the diffusion-coefficient definition
We assess the sensitivity of the diffusive-leakage spectrum to the diffusion-coefficient formulation, and its effect on the performance of the method. For the calculation of the diffusion coefficient, deterministic lattice-physics codes usually rely either on the B 1 methodology (Hebert, 2009) or on one of the conventional transport approximations (Choi et al., 2015): consistent P N , outflow (or outscatter), inflow (or in-scatter). Several works (Choi et al., 2015(Choi et al., , 2017Smith, 2017) have highlighted the weaknesses of approaches other than the inflow transport approximation, that is unanimously credited to have the most rigorous foundation. For Monte Carlo codes, Liu et al. (2016) proved that the optimal method is the already mentioned CMM, equivalent to the inflow transport approximation in the limit of a homogeneous hydrogen slab.
For the benchmark problems investigated in Section 3, we computed the D G distributions with the options featured by the codes SERPENT and APOLLO2-A: the CMM for the former, the B 1 model for the latter, and the outflow transport approximation for both. As an illustrative case, we consider here the assembly with gadolinium rods of Example 4, for which Fig. 29 shows the diffusion coefficient from the CMM and the outflow transport approximation. The leakage spectra determined with the approximate version of the diffusive model are depicted in Fig. 30. In order to focus on the differences caused by the diffusion-coefficient formulation, they have been computed using the reference environmental flux spectra. The spectrum variation predicted by rehomogenization with the diffusive method is also shown.
The deviations between the two approaches have an effect on the leakage especially in the higher-energy region of the fast group (u P 0:7, i.e., E P 0:12 MeV), where the fine-energy diffusion coefficient computed with the CMM is significantly lower than that from the outflow transport approximation. When looking at the spectrum perturbation, the outflow approach is less accurate in the high-energy region. However, it provides a better approximation in the epithermal range, even if the leakage spectra corresponding to the two formulations exhibit negligible differences in this part of the energy domain. In the thermal group, the variations in the spectrum deformation driven by discrepancies in the diffusion coefficient are not significant. If the outflow transport approximation is used instead of the CMM, the errors in R a;1 and R s;1!2 decrease from À0.40% and À0.48% (see calc. e in Table 11) to À0.12% and À0.15% in the Gd-enriched assembly; from 0.39% and 0.45% (see calc. e in Table 10) to 0.1% and 0.07% in the assembly without poison. The error in the total fission power also becomes lower than that found with the CMM formulation: from À1.45% and 1.0% (see calc. e in Table 12) to À0.92% and 0.64% in the fuel bundles without and with gadolinium, respectively. Therefore, the approach used to compute the diffusion coefficient can have a relevant impact on the outcome of the diffusive-leakagebased rehomogenization.  The analysis of the remaining benchmark problems showed that in the thermal group the outflow transport approximation generally provides a more accurate leakage reconstruction. However, as observed for the above sample problem, the impact of diffusion-coefficient variations is small in this energy group. In the fast group, finding an evident trend is difficult and the accuracy of the various formulations depends on the specific test case.

Conclusions and outlook
We presented a method to estimate the interassembly neutron leakage spectrum in the real environment. This computational scheme completes the spectral rehomogenization technique developed at AREVA NP. The proposed approach applies Fick's diffusion law to the node-averaged environmental flux spectra computed by rehomogenization. It uses information from both the fine-group lattice transport calculation and the nodal calculation. The methodology has been extensively validated by numerical simulation of multi-assembly sets exhibiting strong heterogeneity. Several configurations (critical, subcritical and supercritical) have been examined, ranging from fresh-fuel isothermal conditions to depletion feedbacks.
The results show that rehomogenization combined with this leakage model can capture spectral effects on the few-group nodal cross sections very accurately. In the thermal group, the diffusive approach matches the spectral corrections computed with the reference environmental leakage. In the fast group, more accurate cross section inputs for nodal routines are also obtained, even though small deviations from the reference corrections are observed in some configurations. These can be due to: (i) minor inaccuracies in the prediction of the leakage distribution in the high-energy and epithermal regions, and (ii) the difficulty of the basis functions reproducing the resonance details of the spectrum perturbation and its strongly varying shape at high energies. A significant improvement in the estimates of the nodal power distribution and of the multiplication factor is found. The gain in accuracy is achieved at a low computational cost.
The method has been compared to a simpler approach based on the fundamental-mode leakage distribution, computed by the lattice code in the framework of the B 1 critical-spectrum calculation. For all the test cases, the diffusive model outperforms the criticalleakage strategy. The latter only provides a satisfactory approximation of the environmental leakage in the presence of weak spectral interactions, that is, in assemblies surrounded by nodes featuring similar composition. When applied to heterogeneous fuel-loading patterns, this method fails to accurately predict the leakage in the thermal and epithermal ranges of the energy spectrum.
Spectral rehomogenization can only correct a part of the homogenization defect. In order to fully take into account the environmental effects, the method must be applied in combination with spatial rehomogenization. The natural next stage of our investigation is to apply the methodology to reactor problems including thermal-hydraulic feedbacks. Broadening the range of test cases to reflector boundaries is also of great interest. These further steps will be followed by the final validation in full-core calculations. Additional topics to be addressed encompass the impact of spectral rehomogenization on the accuracy of the pin-power distribution reconstructed by dehomogenization and on the calculation of other relevant operational parameters, such as the control-rod bank worth.