Coupling Inward Diffusion and Precipitation Kinetics; the Case of Nitriding Iron-Based Alloys

A model that describes the inward diffusion of an element I into a solid substrate and the simultaneous precipitation of a compound MyIz, with M as the alloying element initially dissolved in the substrate matrix, is presented for the case of nitriding iron-based alloys. The model was developed by coupling the diffusion kinetics and the precipitation (nucleation and growth) kinetics. Additionally, the role of excess nitrogen and the kinetics of ammonia dissociation at the iron surface were incorporated into this coupled model. The model was successfully applied to the case of nitriding an Fe-2.23 at. pct V alloy; the simulation results are in good agreement with the measured data and allow for detailed understanding of the evolution of the nitride precipitates (volume fraction, number density, and size distribution) as a function of both nitriding time and depth in the specimen. The present model exposed the pronounced effects of the precipitation kinetics, of excess nitrogen, and of the surface-reaction kinetics on the overall nitriding kinetics and demonstrated a striking, nonmonotonous change with time of precipitate particle size at a distinct depth in the specimen.


I. INTRODUCTION
THERMOCHEMICAL surface treatments of alloys, such as nitriding, carburizing, nitrocarburizing, carbonitriding, and oxidation, involve uptake and inward diffusion of (interstitial) atoms I (N, C, and O) into the matrix, often in association with simultaneous precipitation of initially present, dissolved alloying element M (Ti, V, Cr, Al, etc.) as nitrides, carbides, and/or oxides (M y I z ). The amount, size, and morphology of the M y I z precipitate particles have a large influence on the mechanical properties of the alloys after the treatment. Therefore, it would be highly beneficial to have at one's disposal a model that links the inward diffusion of I with the (depth-and time-dependent) precipitation state (volume fraction, number density, and size distribution) of the M y I z particles so that the microstructure of the alloys can be predicted and optimized. [1] Until now a series of attempts have been made to simulate inward diffusion of I and simultaneous precipitation of M y I z , [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] usually focusing on one specific (extreme) case of inward diffusion and precipitation. A simple model of coupled diffusion and precipitation was thus developed for the case of internal oxidation [2,3] and later applied to nitriding [4] and carburizing. [5] This model assumes that the solubility of I in the matrix is negligible in the presence of dissolved M. The component I is thus instantaneously consumed by (internal) precipitation of M y I z , which leads to a sharp reaction front between the (completely precipitated) case and the (unreacted, I-free) core and to a reaction rate controlled by the diffusion of I through the M-depleted matrix. The model was generalized by considering a finite solubility of I, leading to some amount of dissolved I ahead of the reaction front during the precipitation; [6] in this variant, the model has been widely applied to carburizing [7][8][9] and nitriding. [10][11][12][13][14][15][16][17][18][19] A serious shortcoming of the type of models introduced earlier [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] is the disregard of the kinetics (the time and temperature dependencies) of the precipitation reaction. Moreover, these models describe the evolution of only the total amount of precipitated M y I z and cannot predict number density and size (distribution) of the precipitates. A relatively more advanced approach for coupling precipitation kinetics and inward diffusion was, for the case of nitriding, proposed in Reference 20 considering nucleation, growth, and coarsening of nitrides by use of a mean radius approach. [21,22] This type of coupled model provides particle number density, mean particle size, and volume fraction of nitrides as a function of depth and nitriding time, but as a result of the mean radius approach, arbitrary assumptions have to be made on the growth/coarsening kinetics. By contrast, particle coarsening arises naturally when the evolution of the size distribution of particles with consideration of the Gibbs-Thomson effect is modeled. [23,24] Although such description of the precipitation kinetics is frequently used for precipitation in homogeneous alloys, to the authors' knowledge, coupling of inward diffusion with such a type of kinetic precipitation model, to describe the evolution of the depth-dependent particle size distribution (PSD), has not been considered until now.
The inward diffusion of I can be strongly affected by its provision via surface absorption: For gaseous nitriding of iron-based alloys, it is known that the kinetics of nitrogen uptake at the iron surface plays an important role in the development of N content-depth profiles. [25,26] In some investigations, this effect was taken into account, [10,12,16,18] but it has been disregarded in many other studies, [11,[13][14][15]17,19,20] where instantaneous establishment of a (equilibrium) dissolved N concentration at the surface was assumed.
Modeling of inward diffusion and precipitation for the case of nitriding iron-based alloys also requires consideration of the phenomenon of ''excess N''; [13][14][15] experiments show that internally nitrided iron-based alloys absorb more N than expected from (1) the precipitation of all alloying elements as nitride of composition prescribed by equilibrium thermodynamics, and (2) the amount of N dissolved in an unstrained ferrite matrix. [27] This additionally absorbed N is called ''excess'' N, [28] which is ascribed to (1) an enhanced solubility of (weakly bonded, mobile) N in the ferrite matrix due to dilation by the elastic accommodation of the nitride-precipitate/ferrite-matrix volume misfit and (2) adsorption of (strongly bonded, immobile) N at the nitride-precipitate/ferrite-matrix interface. [28][29][30] The contribution of these types of differently bonded N to the total N concentration can be obtained quantitatively from N-absorption isotherms. [31] * Yet, most of the diffusion-precipitation coupled models for internal nitriding published until now ignore this distinct role of excess N. [10,12,[16][17][18]20] The purpose of this study is to develop a, in the sense of the earlier discussion, complete model for coupling the inward diffusion of I and the precipitation of M y I z , for the case of nitriding iron-based alloys, which serves as a model for similar other thermochemical surface treatments. The nitride-precipitation kinetics will be described by a time-and temperature-dependent model for nucleation and growth of second-phase particles, based on a description of the evolution of the PSD, as a function of nitriding depth. The model incorporates the roles of the excess N and of the kinetics of ammonia dissociation at the surface. The earlier model by our group [13][14][15] cannot predict the size and number density of precipitates due to ignorance of the precipitation kinetics in the model development. Furthermore, the gradual increase of the surface nitrogen content as a function of nitriding time (which is due to the finite rate of ammonia dissociation at the specimen surface and the finite rate of nitrogen inward diffusion) was not included in the earlier model for the nitriding response. These limitations have now been removed in the newly developed model, which incorporates the time-dependency of the nitrogen absorption at the interface and, in particular, provides information on the evolution of the particle-size distribution of M y N z precipitates as a function of both depth and nitriding time. The developed model is applied to the case of nitriding an Fe-2.23 at. pct V alloy; the calculated results are compared with the available experimental results. A striking, nonmonotonous change with the time of the precipitate particle size at a distinct depth of the specimen is revealed.

II. THEORETICAL BACKGROUND
The inward diffusion-and precipitation-coupled model for the gaseous nitriding of iron-based alloys was developed under the following assumptions and conditions: (1) A thin-sheet substrate geometry is considered. (2) The precipitates are taken as spherical particles. (3) Nucleation of the precipitates complies with classic, homogeneous nucleation theory. (4) Growth of the precipitates is rate-controlled by the diffusion of the alloying element M (the diffusivity of the interstitial element I is assumed to be always much larger). (5) The M y I z -precipitate particles are pure M-I-compounds of fixed stoichiometry (y:z).

A. Inward Diffusion
The inward diffusion of I into the alloy-sheet substrate (one-dimensional diffusion) is given by Fick's second law: where C I is the concentration of I in the alloy matrix at depth x and time t. D I is the intrinsic diffusion coefficient of I in the alloy matrix. Considering the chemical potential gradient as driving force for diffusion, D I can be expressed as: [32] D I ¼D Ã *For the case of the nitriding conditions applied in the current study, although the N concentration at the VN/matrix interface is locally higher than the solubility limit in a pure Fe-matrix, iron nitrides do not precipitate; [14] the N atoms segregated at the VN/matrix interface are likely more stable than the solute N atoms in the matrix.
where D I * is the self-diffusion coefficient of I and / I is the thermodynamic factor of I with the activity a I and the activity coefficient c I of I in the alloy matrix. X I is the mole fraction of component I.** For a dilute ternary system, c I can be described by: [33] ln where c I 0 is the activity coefficient of I for infinite dilution and e I I and e I M are the interaction coefficients for I-I and I-M in the matrix, respectively (the activity coefficient of M, c M , can be expressed analogously.). X M is the mole fraction of M in the matrix.
For solving the diffusion equation, the boundary conditions (BCs) are required. The flux of I provided at the surface of the substrate by the medium surrounding the substrate, J I,s , can (for gaseous nitriding) be described by: [34,35] where b is the rate constant for the reaction by which I atoms are delivered to the surface of the substrate, C I,eq represents the concentration of dissolved I in the solid substrate (at its surface) in equilibrium with the surrounding medium that provides I, and C I,s indicates the actual concentration of dissolved I at the surface of the substrate. The flux J I,s is equal to the flux of I through the surface described by Fick's first law. Thus, the BC at the surface is given by: In view of the symmetry of the sheet substrate, the (net) flux of I through the center plane of the sheet is zero. Thus, the BC at the center is given by: The nucleation rate is taken to comply with classic nucleation theory, [36] and the particle growth rate is assumed to be controlled by long-range diffusion. The nucleation and growth rates are coupled by their dependency on the mean solute concentration in the matrix, as in well-known models for precipitation in a homogeneous alloy. [24] The PSD is represented by discrete particle size classes, each of a certain population (i.e., number) density. In the present case, for class management, a Lagrange-like multiclass approach [22] is used, where the radius evolution of each class is tracked while its population density remains constant. Classes with a particle radius smaller than a certain cut-off radius are removed from the PSD. The evolution of the PSD with time is computed by numerical integration of the nucleation and growth rates: (i) A new class i with radius R i ¢ is created at each time step considering the nucleation in the time interval Dt after time t (Eq. [11]). Its population N i (the number of nuclei per unit volume of the specimen) is given by: For the nucleation rate, dN/dt, see the next subsection.
(ii) The new radii of all other existing classes i at t + Dt, R i (t + Dt), are calculated by: For the growth rate, dR i /dt, see the Section ''II-B-2.'' (iii) The solute concentrations in the matrix are updated after Dt at each time step.

Nucleation
The rate of nucleation of M y I z from dissolved M and I in the alloy matrix is described by adopting classic nucleation theory: [36][37][38][39][40][41] where N is the number of nuclei per unit volume, N 0 is the number of possible nucleation sites per unit volume (i.e., the number of M atoms per unit volume in the matrix), Z is the Zeldovich factor, b * is the absorption rate of solute atoms for a nucleus of critical radius R * , DG * is the activation energy for nucleation (the nucleation barrier), k B is the Boltzmann constant, T is the temperature in Kelvin, and s is the time lag. Considering chemical, interface, and misfit-strain energy contributions in the development of second-phase particles in a matrix DG * at the critical radius R * (=À 2c/(Dg ch + Dg st )) is given by: [1,36] where Dg ch is the chemical Gibbs energy change per unit volume of precipitate during the precipitation and **For the description of thermodynamics and precipitation kinetics, mole fractions will be used as a measure for the amounts of I and M. Mole fractions of I (X I ) and M (X M ) in the alloy matrix can be obtained from the corresponding concentrations (the number of atoms per unit volume) via the relations X I = C I /C T and X M = C M /C T , respectively. C T and C M are the total concentration of matrix atoms and the concentration of M atoms in the alloy matrix, respectively.
For multicomponent precipitates, b * , considering each atomic species i can be given as 4pR Ã2 , [22] where D i is the self-diffusion coefficient of i in the alloy matrix and a is the lattice parameter of the matrix. In the case of the development of M y I z precipitates, b * is approximated by 4pR Ã2 DMXM a 4 recognizing that D M (D I leading to P Dg st and c are the misfit-strain energy per unit volume of precipitate and the precipitate/matrix interface energy per unit area of the interface, respectively. The precipitate particle is assumed to be stable against thermal fluctuations when it has a radius R¢ slightly larger than R * such that DG = DG * À k B T. The value of R¢ corresponding to DG = DG * À k B T can be approximated by: [36] For formation of a stoichiometric phase M y I z from a dilute solid solution of M and I in the matrix, the chemical driving force Dg ch for nucleation can be given by: [33] where a M and a I are the activities of M and I in the alloy matrix, respectively; V M y I z is the molecular volume of M y I z ; and K sp is the solubility product, i.e., K sp ¼ ða The misfit-strain energy Dg st can be given by: [42,43] where e is the linear misfit parameter, l A is the shear modulus of the matrix, and C 6 is given by 3K M y I z =ð3K M y I z þ 4l A Þ, with K M y I z as the bulk modulus of M y I z .

Growth and coarsening
For the case of small supersaturation of the matrix and growth (positive or negative) only limited by diffusion of M in the matrix, the growth rate of M y I z can be expressed by: [44,45] where D M is the diffusion coefficient à of M in the alloy matrix, X M P (=y/(y + z)) is the mole fraction of M in the precipitate, and d ¼ ðy þ zÞV V M y I z À Á is the ratio of the average atomic volumes of the matrix and the precipitate, with Vas the atomic volume of the matrix.
X M int is the mole fraction of M in the matrix at the precipitate/matrix interface and can be regarded as the mole fraction of M in the matrix in local equilibrium with the precipitate of radius R, X eq R M . Analogous to the equilibrium of dissolved M and I with precipitate M y I z in the absence of an interface and misfit strain (Eq. [12]), the corresponding state of local equilibrium recognizing the presence of matrix/precipitate interface (Gibbs-Thompson effect) and of misfit strain is associated with a solubility product K sp R according to: where X eq R I is the mole fraction of I in the matrix and c eq R M and c eq R I are the activity coefficients of M and I in the matrix in equilibrium with the precipitate of radius R, respectively. The solubility product K sp R can be given as (cf. References 46 and 47): This expression can be derived by combination of Eqs. [10] and [12], using the equality of the thermodynamics of nucleation and of growth. [48] Under the aforementioned assumption (4) The precipitates in a class with R > R * (corresponding to X M > X M int ) have a positive growth rate and, thus, increase in size according to Eq. [14]. On the contrary, the precipitates in a class with R < R * (X M < X M int ) have a negative growth rate and, thus, decrease in size according to Eq. [14]. In this way, the coarsening phenomenon is naturally incorporated in the model for precipitate growth.

C. Specific Aspects Related to Excess N Uptake Upon Nitriding Iron-based Alloys
The total N concentration (C N tot ) of internally nitrided Fe-M alloys can be given as C tot is the concentration of N corresponding with precipitation of M y N z , C N str is the concentration of the (mobile [13] ) N dissolved additionally to C N 0 in the strained ferrite matrix as a consequence of elastic accommodation of the precipitate/matrix misfit, [28] and C N int is the concentration of (immobile [13] ) N adsorbed at the precipitate/matrix interface (see Section I).
For the case of spherical M y N z particles, the value of C N int , which depends on particle size, can be given by (cf. Eqs. [7] and [8]): à For the derivation of Eq. [14], the diffusion coefficient D M is assumed to be a constant [44] (compare with Eq. [2]); in the present case, the self-diffusion coefficient pertaining to the nominal concentration of M in the (binary) alloy is employed.
where B is the number of adsorbed N atoms per unit area of the precipitate/matrix interface. The mobile excess N concentration C N str is part of the total concentration of the N dissolved in the strained ferrite matrix; i.e., C N = C N 0 + C N str . The enhanced N concentration in the strained matrix is accounted for by introducing a concentration value C N,eq for the BC at the surface in the inward-diffusion part of the model (see the earlier Section II-A), which is larger than the concentration value C N,Eq 0 of the unstrained matrix. The M y N z precipitate developing in the matrix shows a positive misfit with the matrix (e > 0). In the case of a fully elastic accommodation of the particle/matrix misfit, the increased solubility C N,eq § is given by: [28,29] where V N is the partial atomic volume of N in the matrix and r is the hydrostatic component of the matrix-stress field caused by all existing particles. Introducing the total volume fraction u of the M y N z precipitate particles, r is given by: [28,42,43] To account for the additional misfit introduced by the (particle-size dependent) amount of excess N adsorbed at the particle-matrix interface, Eq. [20a] can be modified to: where misfit parameter e i ¢ and volume fraction u i ¢ of each particle size class i including the N adsorbed at the precipitate/matrix interface are given by: The parameter f describes the extent to which the full misfit due to building out of the lattice of the M y N z particle by the adsorbed N atoms is experienced. [13,14]

III. MODEL IMPLEMENTATION
For the gaseous nitriding of iron-based alloys, as discussed in the previous section, the diffusion model was coupled with the precipitation model incorporating the emergence of excess N. This coupled model was implemented by the following numerical methods and simulation settings.

A. Numerical Methods
To solve the diffusion and precipitation (nucleation + growth/coarsening) equations, the finite difference method (FDM) was employed with the equations discretized in time by using an adaptive time interval Dt (see below) and in space with a depth interval Dx. For solving the diffusion equation by FDM, the Crank-Nicolson method for partial differential equations [49] and the tridiagonal matrix algorithm [50] were adopted. For solving the precipitation equations by FDM, the Euler method for ordinary differential equations [50] was used.
At the jth time step, the concentration of dissolved N in the ferrite matrix (increased, as compared with the j À 1th time step, by inward diffusion of N, excluding the precipitation at the current time step) at every depth step k, (C N ) k j,diff is calculated by Eqs. [1] through [6]. Then the nucleation and growth (including coarsening) of M y N z (i.e., (N i ) k j and (R i ) k j for each class i of M y N z at every depth) are calculated as described in Eqs. [7] through [17]). Classes with (R i ) k j less than the radius corresponding to the size of one M y N z molecular unit are removed from the PSD and the M and N components are added to the solutes content of the matrix. From the calculated (N i ) k j and (R i ) k j , the nitride particle number density N k j , the mean radius R j k , and the volume fraction ðu M y N z Þ j k are calculated as follows: ½26 § The matrix-stress field generated by existing particles leads to a decrease of nitrogen activity in the matrix with respect to the gas phase and, thus, to an increased solubility of nitrogen (Eq. [19]) for constant activity in the gas phase (i.e., constant nitriding potential). This effect is only considered for incorporating mobile excess nitrogen due to the significant effect of the mobile excess nitrogen on the overall nitriding kinetics. [13] .
where V M y N z is the molecular volume of M y N z . The concentrations of dissolved N (C N ) k j and dissolved alloying element (C M ) k j after inward diffusion of N and realization of nitride precipitation are updated by application of mass balances for the N and the alloying element at depth step k as follows: Þ j k , and D(C N int ) k j indicate the changes in concentration of N and the alloying element caused by precipitation of M y N z and the change in N concentration due to N adsorption at the precipitate/matrix interface during the jth time step at the kth depth step, respectively. The value of C N,eq for the surface BC in the N diffusion model (Eq. [5]) is updated by Eqs. [19] through [22]. The total N concentration at depth step k, (C N tot ) k j , is calculated by: At the (j+1)th time step, N diffusion and nitride precipitation calculations are performed with the parameters updated at the jth time step, and so on. The discussed scheme of calculating N diffusion and nitride precipitation is illustrated in Figure 1.
The shorter the time interval Dt, the better the accuracy of the calculation. However, a too small Dt increases the calculation time without gain in accuracy. Therefore, Dt is adapted during the sequential calculations as follows: Dt is increased by factor 1.1 (Dt fi 1.1Dt) if the critical radius R * updated at the current time step differs less than 1 pct from the value updated at the previous time step for all depths where the precipitation occurred. If this condition is not satisfied, the calculation for the current time step is repeated with a smaller time interval (Dt fi 0.5Dt) with the parameters updated at the previous time step.

B. Simulation Settings
The inward diffusion and precipitation-coupled model developed in this study was applied to the gaseous nitriding of Fe-2.23 at. pct V alloy specimens, with a thickness of 1 mm, at 793 K, 823 K, and 853 K (520°C, 550°C, and 580°C) for 10 hours by using a nitriding potential (r N ) of 0.103 atm À1/2 , i.e., the same conditions as applied in the experiments of Reference 14, so that experiment and simulation results can be compared. Under these nitriding conditions, no iron nitrides can be formed [14] and vanadium nitrides are precipitated in the form of VN. [14,31] This was shown for several Fe-V alloys, where after so-called de-nitriding of specimens previously nitrided until saturation (i.e., no further nitrogen uptake occurs), the amount of remaining, strongly bonded N is at the level of the amount of V in the specimen, indicating that the precipitate composition can be described as VN. [31] Gaseous nitriding of iron-based alloys in a flowing NH 3 /H 2 gas mixture allows precise control of the applied chemical potential of N by the so-called nitriding potential r N (=p NH 3 . p 3=2 H 2 , p is the partial pressure) [51] and of the kinetics of the dissociation of ammonia at the iron surface. [27,34,35,52] The diffusion and precipitation-coupled model including the effect of excess N and surface reaction kinetics is termed full model. Simulations were also performed with the model in a less ''complete'' state, to expose the role of various, individual processes integrated in the full model: (1) To illustrate the importance of considering the role of excess N, simulations neglecting the role of excess N have been performed; the model in that fashion is called the basic model. (2) To investigate the effect of the time (and temperature) dependence of the nitride-precipitation kinetics, simulations were also performed under the assumption of instantaneous establishment of local thermodynamic equilibrium in the diffusion zone, implying instantaneous nitride precipitation once the solubility is (locally) surpassed. This is called the ''equilibrium approach'' (EA) (see Section I). (3) To illustrate the effect of neglecting the kinetics of N absorption at the surface, simulations were carried out assuming instantaneous establishment of local equilibrium with the NH 3 /H 2 gas mixture at the surface of specimen; i.e., the following BC was used: The simulations were performed with a depth interval Dx of 10 lm (and with a smaller Dx of 1 lm especially for the time evolutions of the precipitation state according to the basic model and for the EA model, which both are numerically unstable for large Dx) and an initial time interval Dt of 0.01 seconds.
Thermodynamic data for c N , c V , and K sp for VN, which were used in the present simulations, were obtained by using Thermo-Calc software and the TCFE7 database, [53] and they have been listed in Table I. Other used input data are presented in Table II. The temperature dependence of the bulk modulus of VN was assumed to be identical to that of TiN for which literature data are available. [62] The value of B (cf. Eq. [18]) was obtained by using the measured value of C N int (0.40 at. pct) for homogeneously, through nitrided Fe-2.23 at. pct V thin foil (thickness of 0.2 mm) specimen, nitrided for 24 hours at 853 K (580°C) at r N = 0.103 atm À1/2 . [31] The value of B was taken as temperature-independent.

IV. RESULTS AND DISCUSSION
A. Basic Model (i.e., Ignoring Excess N) Results of the basic model, for the case of nitriding Fe-2 at. pct V alloy at 853 K (580°C) for 10 hours using  a nitriding potential of 0.1 atm À1/2 , are shown in Figure 2 in terms of the contents of dissolved N and V in the matrix, the content of N present in the form of VN-particles, and the volume fraction, the number density, and the mean radius of the VN-particles, all as a function of depth below the surface. For this nitriding time, up to a depth of 250 lm, nearly all V has precipitated as VN (Figure 2(a)), which is expected due to the negligible solubility of V in the ferrite matrix at the employed nitriding conditions. The small negative depth gradient of the total N content in the diffusion zone (at x <~250 lm) originates from the depth gradient of N dissolved in the matrix. In accordance with a constant solubility product of nitride, the content of V dissolved in the matrix shows a positive depth gradient (see the inset of Figure 2(a)); i.e., the dissolved V content increases with decreasing dissolved N content. Within the region of (practically) completed VN precipitation (x <~250 lm), an interesting variation of particle number density and mean particle size as a function of depth occurs (Figure 2(c)): The number density of VN particles is very high at the surface, but it decreases rapidly with increasing depth (note the logarithmic scale in Figure 2(c)). In accordance with the practically constant volume fraction of VN within this region (Figure 2(b)), the mean radius of the VN particles shows an opposite trend. These observations indicate a strong variation of the rates of nucleation and growth as a function of depth. Due to the continued, fast provision of N by the decomposition of ammonia at the surface, the region adjacent to the surface of the specimen is exposed to a higher flux of N than at distinct depths, where N has to be provided by solid state diffusion of N from the surface. The inward flux is additionally reduced by the precipitation of VN in regions closer to the surface. Consequently, the driving force for nucleation in the region adjacent to the surface is higher than in the subsurface regions, leading to a higher number density and smaller size of VN particles in the surface-adjacent regions. This phenomenon was indeed experimentally observed in transmission electron microscopy (TEM) investigations carried out on nitrided Fe-Cr [63] and Fe-Cr-Al [64] alloys.
The earlier discussed effects are even more dramatically illustrated considering the time evolutions of the number density, size, and volume fraction of VN at the surface (depth = 0 lm) and at the depth of 200 lm below the surface (Figure 3).
At the surface, the number density of VN particles increases rapidly after the onset of nitriding (Figure 3(a)). The volume fraction of VN, however, is still low (Figure 3(c)) due to the very small size of the VN particles (Figure 3(b)). With increasing nitriding time, the number density of VN particles remains almost constant (after about 6 seconds), but a gradual increase in mean radius and amount of V precipitated as VN occurs due to continued particle growth (almost) up to the equilibrium fraction (at about 900 seconds). Hence, two different stages of the precipitation reaction at the surface can be distinguished: an initial stage dominated by rapid nucleation (stage I), followed by a stage of particle growth with negligible nucleation, i.e., almost pure growth (stage II). This behavior is caused by the high sensitivity of the nucleation rate on the degree of supersaturation, which is determined by the contents of dissolved N (given by competing rates of its provision at the surface and its consumption by the precipitation) and dissolved V in the matrix. In stage II, a low level of supersaturation remains due to the relatively high consumption rate of N by growth of the large number of existing particles. With increasing time, in the surface-adjacent region, no significant change in the number density of particles and the mean radius is observed indicating negligibly slow coarsening kinetics (stage III): The coarsening rate is very low due to the extreme depletion of V in the matrix and the continuous supply of N, which, in view of the solubility product of VN, means that the rates of both positive growth of larger particles and negative growth of smaller particles are very slow (Eqs. [14] and [17]). The evolution of the PSD of VN at the surface is illustrated in Figure 4; after 1.5 seconds (in stage I), after 20 seconds (in stage II), after 2000 seconds (in stage III), and after 36,000 seconds, i.e., after the total nitriding process time of 10 hours. The results shown reflect the earlier interpretation of the evolution of the precipitation state in the surface-adjacent region: The particle number density increases in stage I (Figures 4(a) and (b)); the particle size increases in stages I and II (Figures 4(a) through (c)); and then the PSD broadens in stage III (Figures  4(c) and (d)).
At the depth of 200 lm, the evolution of the state of precipitation can also be subdivided into three stages, as shown in Figure 3. Stages II and III are similar to those described for precipitation in the surface-adjacent regions, i.e., almost pure growth in stage II and slow  coarsening in stage III. However, stage I differs strongly from that observed for the surface-adjacent region: the mean radius initially increases continuously with time, but this is then followed by a decrease toward the end of the stage I, in contrast with the continuous increase in mean radius occurring at the surface until the end of stage I. This phenomenon can be explained as follows: Due to the much lower flux of N as compared with the surface-adjacent region, the degree of supersaturation attained at the moment of VN precipitation and, thus,  the driving force for VN precipitation is pronouncedly smaller at a distinct depth in the diffusion zone than in the surface-adjacent region. Thereby a much smaller number of particles (Figure 3(a)) of larger size nucleate and then grow (Figure 3(b)) in the early stage of the VN precipitation, in association with a relatively low consumption rate of N as compared with the surface-adjacent region. With increasing nitriding time, upon (virtual) completion of the VN precipitation in the surface-adjacent region, at a depth of 200 lm, the provision rate of N by inward diffusion gets larger than its consumption rate by VN precipitation, leading to an increase of the dissolved N content in the matrix; the driving force for nucleation at this depth can then increase despite the already reduced dissolved V content in the matrix. This results in the nucleation of a relatively large number of particles (Figure 3(a)) of size smaller than the momentary mean radius, leading to a decrease of the mean radius (Figure 3(b)). After this more or less abrupt increase of the nucleation rate, at a next moment of time at this depth, the nucleation rate decreases due to a decrease of the driving force for nucleation by the consumption of the solutes V and N, and particle growth becomes relatively dominant causing a visible increase of mean radius: Stage II of almost pure growth is entered.

B. Comparison with the Equilibrium Approach Model
According to the EA model, at any location all N above the equilibrium content corresponding to the solubility product of the nitride is consumed for instantaneous precipitation of nitrides. As follows from the results shown in Figure 5(a), the EA approach, as expected, leads to a sharp transition between the nitrided case and the unreacted core, whereas incorporating nitride-precipitation kinetics (the basic model) leads to a more gradual decrease in the total N content with a depth between the fully precipitated case and the unreacted core. As shown in Figure 5(b), in the basic model, the overall evolution rate of the amount of VN is much lower. This allows a larger inward diffusion flux of N into deeper regions, which results in an earlier onset of precipitation in these deeper regions than for the EA model.

C. Comparison with the Model Ignoring Surface Reaction Kinetics
Disregarding the ammonia dissociation kinetics at the surface implies. C N | x=0 = C N,eq for t > 0. Consequently, a larger inward flux of dissolved N occurs (Eq. [5]), especially in the early stage of nitriding. This leads to larger values of C N tot at every depth and, in particular, to a more pronounced nitriding depth as shown in Figure 6. In reality, as indicated by the results obtained from the basic model, C N at the surface only equals C N,eq after through thickness, homogeneous nitriding of the specimen. D. Full Model (i.e., Incorporating Excess N) The N content-depth profiles of Fe-2.23 at. pct V alloy nitrided with r N = 0.103 atm À1/2 at three different temperatures [793 K, 823 K, and 853 K (520°C, 550°C, and 580°C)) for 10 hours, as measured [14] and as simulated according to the full model and the basic model, are shown in Figure 7. Evidently, incorporating the role of excess N (mobile and immobile) is essential to arrive at results compatible with experimental observations: The results obtained by the basic model show significant disagreement with the measured data, whereas the results obtained by the full model match fairly well with the experimental data at all temperatures.
Accounting for the effect of excess N greatly influences the predicted nitriding kinetics: The mobile (=dissolved in the matrix) excess N, C N str , mainly induces an increase of the nitriding depth and slightly contributes to an increase of the total N content, C N tot , at every depth where nitride precipitation occurs [13] by enhancing the flux of N diffusion into the bulk. (Eqs. [5] and [19] through [22]). The immobile (=adsorbed at the precipitate/matrix interface) excess N, C N int , primarily causes an increase of C N tot at every depth where nitride precipitation occurs, and it induces a small decrease of nitriding depth [13] by consumption of N in the region of precipitation.
The results of the full model and the basic model for the nitriding at 853 K (580°C) in terms of VN number density, VN mean radius, VN volume fraction, V content dissolved in the matrix, N content dissolved in the matrix (C N = C N 0 + C N str ), N content consumed for precipitation of VN, and N content adsorbed at the precipitate/matrix interface (C N int ) are shown as a function of depth in Figure 8. The trends as a function of depth according to the full model as shown in Figures 8(a) through (f) are similar to those observed for the basic model (see our earlier discussion). However, the quantitative results are different, as a consequence of the incorporation of the role of excess N: Compared with the basic model, the full model predicts larger particle density and smaller particle size. This is due to the larger driving force for nucleation caused by the larger N flux from the surface into the bulk, which is induced by a higher value of C N,eq (Eqs. [5] and [19] through [22]) due to the misfit-strain field incorporated in the full model. The particle size calculated with the full model (the mean radius, in the region from the surface to the depth of 200 lm, equals 1.2 to 1.9 nm; the corresponding average particle volume varies from 7.2 to 28.7 nm 3 ) is compatible with the very tiny VN particles (5 nm long and one-to two-monolayer-thick platelets, with a corresponding average particle volume of 5.2 to 10.4 nm 3 ) as observed in TEM investigations performed on the same alloy specimens of different thicknesses (200 lm) nitrided under identical nitriding conditions. [65] It is of interest to note that C N int decreases with increasing depth in the precipitated zone (Figure 8(g)) despite a practically constant VN fraction in this region (0 to 380 lm), where nearly all V has precipitated as VN. This is a consequence of the decrease of precipitate/matrix interface area (per unit volume of VN) due to the increase in mean particle radius with increasing depth.
V. CONCLUSIONS 1. A general model for describing the coupling of the inward diffusion of an element I and the simultaneous precipitation of a compound M y I z , with M as the alloying element initially dissolved in the matrix, was developed. The diffusion submodel incorporates the competition of the surface reaction, providing I to the surface of the substrate, with the inward diffusion of I. The precipitation submodel details the concurring processes for nucleation, growth, and coarsening, incorporating the Gibbs-Thomson effect and the effect of misfit strain. The coupled model allows for predicting the amounts of I and M in the matrix and in the precipitate particles, and the state of precipitation, as given by particle volume fraction, number density, particle mean radius, and particle size distribution, all as a function of time, temperature, and depth. The coupled model can be applied to cases of internal nitriding, carburizing, oxidizing, etc. 2. The coupled model was applied to the gaseous nitriding of an Fe-2.23 at. pct V alloy. The calculated N content-depth profiles and particle sizes show good agreement with experimentally obtained data. 3. The simulation results, obtained for the case of nitriding, provide detailed understanding of the coupling of the inward diffusion and precipitation reaction: (a) With increasing depth, the particle number density decreases and the mean particle radius increases due to the decrease in driving force for nucleation caused by the decreasing flux of the inward-diffusing element, here N. (b) Depending on the degree of supersaturation, either nucleation (high supersaturation, stage I), growth (lower supersaturation, stage II), or coarsening (very low supersaturation, stage III) can become dominant.
(c) For regions remote from the surface, the intensity of the flux of the inwardly diffusing element is additionally influenced by the state of precipitation in regions closer to the surface: As precipitation is close to completion in the surface-near regions, deeper regions can experience an increase of the mentioned flux, leading to an increase in the driving force for nucleation. Thus, at these depths, an enhancement of the nucleation rate of particles occurs; the then nucleating particles have a smaller size than the already present ones, leading to a decrease of the mean particle radius. 4. The separate effects of the nitride precipitation kinetics, the surface reaction kinetics, and the excess N uptake are revealed by ''switching on/off'' of these processes in the simulations: (a) Nitride precipitation kinetics leads to a more gradual case/core interface than instantaneous precipitation upon surpassing the solubility limit. (b) The surface reaction kinetics significantly reduces the overall nitriding rate in association with the eventual establishment of a stationary state at the surface.
(c) The presence of excess N greatly enhances the overall nitriding kinetics and total N content at each depth; it induces an increase in nitride-particle number density and a decrease in nitride-particle size.
ACKNOWLEDGMENT Open access funding provided by Max Planck Society (Institute for Intelligent Systems).

OPEN ACCESS
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.