Modelling the reduction of quartz in a quartz–carbon pellet

Traditional refining of silicon generates carbon dioxide emissions that are released into the atmosphere. An alternative process under experimental exploration uses quartz particles coated in a layer of porous carbon as the raw material, instead of lumps of quartz and carbon. The quartz–carbon pellets are processed at a lower temperature than in an industrial furnace, and hence, different chemical reactions are dominant, reducing the greenhouse gas emissions. The quartz core shrinks as it is consumed, and the carbon is converted to silicon carbide, which can subsequently be processed into silicon. We develop a model for chemical and transfer processes within a single quartz–carbon pellet. We derive governing equations for the concentration of silicon monoxide, carbon monoxide, and carbon dioxide, and conservation equations on the moving quartz interface. Furthermore, we then focus on a reduced model in an industrially relevant distinguished limit, and solve numerically the resulting leading-order system. We show examples of reaction-limited behaviour as well as diffusion-limited behaviour. Both regimes are physically admissible due to the large potential range of the parameters. Finally, we sweep through the parameter space, and characterize the dynamics based on the utilization of the carbon and the silicon yield. We find that the diffusion-limited regime is best for carbon utilization and silicon yield, as the silicon monoxide reacts with the carbon before it is transported out of the pellet.


Introduction
The second most abundant element in Earth's crust is silicon [1], occurring naturally as silicon dioxide in quartz. Typically, silicon is obtained by processing quartz, or sometimes quartzite [2]. Silicon is used as a semiconductor for use in electronics such as solar panels, and is also used to produce silicone and alloys [2]. Quartz is often processed into silicon in a submerged-arc furnace, and similar furnaces are also used for producing ferrosilicon, phosphorus, and magnesium [2].
Carbon monoxide is the primary by-product of the traditional silicon process in a submerged-arc furnace. The carbon monoxide combusts to carbon dioxide, and is normally released into the atmosphere [3]. Metal and metalloid factories are among the most intensive and concentrated sources of greenhouse gas emissions in the world [4]. Typically, about 6.1−6.5 kg of carbon dioxide is released into the atmosphere for every 1 kg of silicon produced [2]. Since emissions involved in the production of the carbon sources that are used in the furnace is not accounted for, the figure is an underestimate. Alternate processes are therefore being explored to reduce, or neutralize, the carbon footprint of silicon production.
One alternative process being explored involves granular quartz-carbon pellets, shown in Fig. 1, as a precursor. The carbon in the pellets is converted to silicon carbide, which can then later be processed into silicon. Using quartz-carbon pellets potentially improves the energy efficiency, and may significantly decrease greenhouse gas emissions [4]. The pellets are created by depositing carbon from renewable methane onto a quartz core [5]. Since the carbon layer is deposited by a chemical reaction, the carbon is extremely pure [4,5], unlike more traditional forms of carbon used in furnaces, and higher purity carbon yields higher purity in the final silicon [5]. The deposited carbon forms a porous layer surrounding the solid quartz core [6]. We show an electron microscope image of the carbon layer in Fig. 2 to highlight the porous structure. In a typical laboratory-scale experiment, the pellets are added to a cylindrical crucible. The pellets are then heated up, starting from room temperature, by an inductive furnace [7]. The temperature is then held at around 1600 • C for a specified amount of time, about an hour [6], and then the inductive heating is turned off and the crucible is allowed to cool down [6]. By using lower temperatures than usually exploited in traditional furnaces, different chemical reactions can be made to be dominant. After the experiment, the crucible is injected with epoxy. Cutting open the solidified slab of epoxy reveals the exact cross-section of the crucible at the end of the experiment. The slab can be visually inspected for consumption of the reactants and formation of the products and differences throughout the crucible. Additionally, samples can also be analysed with an electron microscope or with X-ray spectroscopy [6][7][8] to understand the microstructure.
We briefly review the traditional silicon process [2] before returning to the chemistry relevant to quartz-carbon pellets. A traditional silicon furnace is continually filled with a mixture of solid quartz and a form of carbon-usually a combination of charcoal, coal, woodchips, and coke. The mixture of quartz and carbon is referred to as the charge. A large, three-phase alternating current passes between three electrodes and the bottom of the furnace. The current heats the quartz and carbon owing to Ohmic resistance, and the quartz begins to melt, initiating the chemical reactions. The main reactions within an industrial furnace, once the quartz has melted, are [9] SiO 2 (l) + C(s) − − → SiO(g) + CO(g), SiO(g) + SiC(s) − − → 2Si(l) + CO(g).
Initially, the quartz (SiO 2 (l)) and the carbon react to form silicon monoxide (SiO) and carbon monoxide (CO). The silicon monoxide then reacts with carbon to form silicon carbide (SiC) and additional carbon monoxide, or condense to form liquid quartz and silicon (Si(l)). Finally, the silicon monoxide can also react with silicon carbide to yield liquid silicon and carbon monoxide gas. The liquid silicon is tapped off at the side of the furnace, while the gaseous silicon monoxide and carbon monoxide are vented through the charge. As the gases rise through the charge, they begin to cool; some silicon monoxide condenses on the surface of the charge and forms a viscous liquid of quartz and silicon via reaction (3). The remaining silicon monoxide and carbon monoxide combust once outside the furnace into silicon dioxide and carbon dioxide, respectively. As the silicon dioxide and carbon dioxide cool further, the silicon dioxide condenses into microsilica, which is captured and sold as a concrete strengthener [2,3], whereas the carbon dioxide is released into the atmosphere [2]. The kinetics and chemistry within a standard furnace have been widely studied, for example [9][10][11]. Similarly, particular sections of a furnace have been modelled mathematically, including the arc and current distribution [12], the evolution of the cavity and the charge [7], and the formation of microsilica as silicon monoxide combusts [3]. We now consider the quartz-carbon pellets used in the experimental exploration. Here, the quartz is below its melting temperature of ∼ 1700 • C, and the chemical reactions are different from (1)- (4). The initiating reaction is uncertain [9], although Li [6] suggests is the initial reaction. We notice that reaction (5) is almost the same as reaction (1). Reaction (5) can happen only for a short time since it is a solid-solid reaction-once the first layer of atoms have reacted, the two solids lose contact. The silicon monoxide formed via reaction (5) then reacts with the carbon layer, by reaction (2), as the silicon monoxide is transported out of the pellet. Conversely, carbon monoxide reacts with the surface of the solid quartz core via The final reaction for quartz-carbon pellets is the Boudouard reaction, which may be limited by the production of carbon dioxide via reaction (6). Reactions (2), (6), and (7) may be combined to yield the net reaction Reaction (8) has been widely studied in the literature [10,[13][14][15][16][17][18][19][20][21][22][23]. The ideal carbon to quartz mole ratio is 3. However, to maximize the silicon yield, we would want to fully consume the initial quartz. Thus, it is common experimentally to use an excess of carbon-a mole ratio of 3.2-3.5 typically. The reactions (2), (6), and (7) were initially proposed as the mechanism for the reduction of quartz using carbon by Lee, Miller, and Cutler [19]. In our mathematical model we will include these reactions, the silicon Fig. 3 Cross-section of a single, spherically symmetric, quartz-carbon pellet. A quartz core is initially surrounded by a porous carbon layer and the gaseous exterior. The quartz core shrinks owing to reaction (6), creating an internal gaseous region in the pellet between the quartz and the carbon monoxide-carbon reaction, (2), the carbon monoxide-quartz surface reaction, (6), and the Boudouard reaction, (7). The quartz-carbon pellets change in shape and size during heating, and has been studied by Li [6] and Li et al. [24]. Similar composites have also been studied [6,13,[25][26][27][28]. These composites are created by binding ground (∼ µm diameter) quartz and carbon into spherical particles (∼mm-cm diameter). The key difference between the pellets and composites is that the composites are a homogeneous mixture of quartz and carbon, whereas the quartz-carbon pellets we are interested in have the discrete layers. Understanding the production of silicon carbide and the escape of silicon monoxide in composites has been well studied [6,13,27,28]. Generally, the conversion of the quartz and carbon to silicon carbide depends on the mole ratio of quartz and carbon, the size of the ground particles and the composites, and the pressure used when creating the composites. Finally, Li and Tangstad [25,26] have compared the dependence of silicon carbide formation on temperature, carbon monoxide gas concentration, and carbon content for the pellets and the composites.
In this paper, we will model the mass transfer within a single quartz-carbon pellet to estimate the consumption of the carbon and the conversion of quartz to silicon carbide to predict the efficiency. We will present a mathematical model of the chemical reactions and transport, and then non-dimensionalize the model in Sect. 2. We will then reduce the model by considering a physically relevant regime in Sect. 3. The resulting simplified model allows us to partially solve the model analytically. We will explore the simplified model to understand the reaction-limited and diffusion-limited behaviour in Sect. 4. We will then investigate sweeps of the parameter space to better understand the transition between the reaction-limited and diffusion-limited behaviour. Furthermore, we will quantify the behaviour by evaluating how much of the initial carbon is consumed, as well as how much of the initial quartz is converted to silicon carbide. Finally, we will summarize our results and give our concluding thoughts in Sect. 5.

Mathematical model
We consider an idealized case of an isolated spherical quartz particle with a concentric, porous carbon shell surrounding the quartz, as shown in Fig. 3. Our geometry consists of three spatial domains. Initially, a solid, impermeable quartz core of radius r = R 1 is surrounded by a porous carbon layer, extending from r = R 1 to r = R 2 . We neglect the effect of gravity, since we wish to maintain the critical spherical symmetry. We assume all regions initially contain an inert gas (such as nitrogen or argon [9,10,22,28]) and some externally supplied carbon monoxide with a concentration of c 0 . In order to sustain the reactions, we continue to provide the external source of carbon monoxide in the far-field. As time progresses, the carbon monoxide-quartz reaction, (6), consumes the quartz core-reducing the size-and we denote the resulting radius of the quartz core by r = s(t). As the radius of the quartz shrinks, an interior gaseous cavity is formed in the region s(t) < r < R 1 .
We assume the gas mixture contains four components: silicon monoxide, carbon monoxide, carbon dioxide, and inert nitrogen. We denote the gaseous concentration of silicon monoxide by C SiO , carbon monoxide by C CO , and carbon dioxide by C CO 2 , each with units of mol m −3 of gas. Furthermore, we assume that the gas is dominated by the inert gas, and so we assume the gas mixture has constant density.
For simplicity we assume the system is isothermal, since the pellet's relatively small size suggests that the temperature will be uniform throughout the pellet, and for simplicity we assume that the temperature remains constant.
In the exterior of the pellet, r > R 2 , the assumed spherical symmetry implies conservation of mass for each of the reactive species reads where the fluxes q E i are defined to be for i ∈ {SiO, CO, CO 2 }; here u E is the radial mass-averaged speed of the gas mixture, and D i is the diffusivity of species i in nitrogen, which we take to be constant. The superscript E on the dependent variables is to make it clear these quantities are defined in the exterior region, r > R 2 . Instead of directly modelling the concentration of the inert gas, we appeal to conservation of total mass of the gas and the assumption of constant density, and write down In the gap between the quartz and carbon, s(t) < r < R 1 , which we call the gas cavity, we have the same dynamics as in the exterior, and hence, our model in the gas cavity is identical to (9)-(13), but with the superscript I , for interior gas cavity, labelling the dependent variables.
Within the porous layer the silicon monoxide-carbon reaction, (2), and the Boudouard reaction, (7), take place. The expressions representing conservation of mass in the porous layer, R 1 < r < R 2 , are given by where φ is the porosity of the porous layer, ρ 0 is the density of the gas, and the intrinsic fluxes q P i are defined to be The superscript P denotes these quantities are defined in the porous layer. Since C P SiO is measured in mol m −3 of gas, φC P SiO is the extrinsic concentration of silicon monoxide per unit volume of porous media in mol m −3 . Here, R 1 and R 2 are the reaction rates of the silicon monoxide-carbon reaction, (2), and of the Boudouard reaction, (7), respectively-we will specify constitutive laws for the rates in Sect. 2.2-and D p i is the diffusivity of species i in the void space of the porous medium, which we again assume to be constant. On the scale of the pores within the porous medium, we consider reactions (2) and (7) to be surface reactions. However, since the reactions take place throughout the porous medium, we consider R 1 and R 2 to be bulk reactions on the scale of the porous medium.
Within the solid phase of the porous layer, we assume the carbon and silicon carbide are immobile. We assume the thickness of the porous media is constant, and any volume change due to reaction (2) results in a change in porosity instead of an increase in the thickness of the porous medium. It then follows that conservation of mass simply reads where C C and C SiC are the concentration of carbon and silicon carbide in mol m −3 of solid, respectively. We note that (1 − φ)C C is the extrinsic concentration of carbon per volume of porous media and since the solid is made of only carbon and silicon carbide, we have where m C and m SiC are the molar mass of carbon and silicon carbide, respectively, and ρ C and ρ SiC are the density of carbon and silicon carbide, respectively. Multiplying (21) by 1 − φ, and taking the time derivative, we find by (19) and (20). The R 1 term represents the difference in volume caused when two carbon atoms are replaced by one silicon carbide molecule, while the R 2 term is the volume lost owing to the Boudouard reaction. We neglect any structural issues that may arise in the limit φ → 1. There may be some maximum porosity, φ max < 1, above which, the carbon layer disintegrates. In our numerical implementation, we will use (22) in place of (20), and use (21) to find the concentration of silicon carbide a posteriori.
We note that we assume the thickness of the porous layer is constant. As the carbon is consumed, some of the carbon is replaced by silicon carbide. As a result, the region R 1 ≤ r ≤ R 2 maintains the same thickness, with a different porosity. Thus, this region is porous with φ < 1, even if the carbon consumption is concentrated at the boundaries. Initial experimental evidence has shown the thickness does not vary.

Boundary conditions
We recall that we denote the boundary of the solid quartz core by r = s(t), and the inner and outer boundary of the porous layer by r = R 1 and r = R 2 , respectively. As r → ∞, we assume the concentration of carbon monoxide approaches a constant, and the concentrations of silicon monoxide and carbon dioxide tend to zero, so that where c 0 is the concentration of the external supply of carbon monoxide. On the interface between the exterior and the porous layer we assume the concentrations are continuous, so that We also assume continuity of flux of silicon monoxide, carbon monoxide, carbon dioxide, and mass, so that Similarly, on r = R 1 , the boundary between the gas cavity of the pellet and the porous layer, we have analogous boundary conditions to (24) and (25).
On r = s(t), the quartz interface, we consider the effect of the carbon monoxidequartz reaction, (6). We assume no gas enters the solid quartz core, and so the five boundary conditions needed are obtained from considering conservation of the three gas species, the total mass of the gas, and the mass of the quartz. Conservation of silicon monoxide relative to the moving interface yields where R 3 is the reaction rate of the carbon monoxide-quartz reaction, (6), to be specified in the following subsection. Similarly, we have the boundary conditions for conservation of carbon monoxide, carbon dioxide, and the total gas mass, respectively, where ρ SiO 2 is the density of quartz. We note that the right-hand side of (27) is negative, since carbon monoxide is consumed by reaction (6). The rate that the quartz is consumed is related to the carbon monoxide-quartz reaction, (6), by where m SiO 2 is the molar mass of quartz. Unlike the constant thickness of the porous layer, the quartz decreases in size. By producing solid silicon carbide, the thickness of the porous layer remains constant, whereas the products of the (6) are all gaseous, and therefore the quartz shrinks.

Rates of reaction
We now specify the form of the rates of reaction R 1 , R 2 , and R 3 . The particular rates of reaction involving quartz and carbon are complicated and not well understood [9]. For simplicity, we assume reactions (2), (6), and (7) are elementary reactions, and model the reaction rates with the law of mass action by assuming where k 1 and k 2 are the associated rate constants, with units m 3 / (mol s). In reality, k 1 and k 2 will be temperature dependent, however, we assume they are constant by our isothermal assumption. We note that we assume R 1 and R 2 are linear in C C since the rate limiting step is the gas reaching the solid, see [29]. The carbon atoms are densely packed compared to the gas, and so, whenever a silicon monoxide molecule interacts with a carbon atom, there is always an adjacent carbon atom to allow the silicon monoxide-carbon reaction, (2), to proceed. For the carbon monoxide-quartz reaction, (6), we assume the form where k 3 is the associated rate constant, with units m / s, and H (x) is the Heaviside function. We do not explicitly track the quartz concentration in (33) since it is constant, and instead is contained in k 3 , as is usually done. Once again, in reality k 3 will depend on the temperature, but we assume k 3 is constant. Note, the carbon monoxide-quartz reaction, (6), is a surface reaction, whereas reactions (2) and (7) are bulk reactions. Hence, the Heaviside function is required to turn off the reaction once the quartz has been fully consumed, and we transition from modelling a quartz-carbon pellet to modelling a hollow carbon shell. Once we are in the latter regime, the partial differential equations that were imposed for s(t) < r < R 1 must now be imposed on 0 < r < R 1 . The boundary conditions on r = s(t) are then replaced in our symmetric geometry by the condition that the concentrations are bounded as r → 0 + , and that u I = 0.

Initial conditions
Initially, the quartz core extends to the porous carbon layer, and so The solid is taken to be composed entirely of carbon. Thus, at t = 0 we have We assume no silicon monoxide or carbon dioxide is present in the system initially and the carbon monoxide concentration is uniform, and thus write at t = 0; we note that these initial concentrations are consistent with the far-field boundary conditions given in (23). We now consider typical sizes for our various parameters; the values are given in Table 1. The uncertainty of the value of each parameter varies considerably from parameter to parameter. Physical constants, such as the density of carbon or the molar mass of silicon carbide, are accurately known. However, the rate constants are not known experimentally, and their values must be estimated from related experiments. We expect the physical range of the rate constants to be much narrower than reported   Table 1, however, we consider this wide range to ensure we capture the correct value in our simulations. In the experiments, the pellets have a total diameter of about 400 µm. However, the radius of the quartz core and the thickness of the carbon layer can be varied. We wish to explore the dynamics over a range of R 1 and R 2 values. The initial carbon monoxide concentration can also easily be varied; we impose a 10% molar concentration.

Non-dimensionalization
We non-dimensionalize the model (9)-(36) in the following way. We scale r and s with the outer radius of the pellet, In each region we scale C SiO , C CO , and C CO 2 by the initial concentration of carbon monoxide, and write Multiple possible time-scales are embedded in our model. We choose to scale time with the conversion of carbon to silicon carbide, found by balancing terms in (19). We choose the velocity scale to balance the flow generated on the surface of the quartz with the speed of the quartz interface in (29), noting that ρ SiO 2 ρ 0 . We thus write Finally, we scale the concentrations of carbon and silicon carbide by the initial concentration of carbon, and write After dropping the tildes, the dimensionless model in the exterior, from (9)-(13), is for r > 1, where Here, q E i are the dimensionless fluxes, and D i and Pe i are the dimensionless diffusivities and Péclet numbers.
In the interior gas cavity we have analogous equations to (41)-(44) which hold in the dimensionless region s(t) < r < s 0 , where In the porous layer, the dimensionless equations for the concentration of carbon and silicon carbide, and the porosity, (19), (20), and (22), respectively, are which hold in the dimensionless region s 0 ≤ r ≤ 1. We have introduced the two dimensionless parameters which are the relative volume change when two carbon atoms are replaced by one silicon carbide molecule, and the ratio of the rate constants for the silicon monoxidecarbon reaction, (2), and the Boudouard reaction, (7), respectively. The equations governing the gas concentrations, (14)- (16), and continuity of mass, (17), become where We have introduced five additional dimensionless parameters. Firstly, δ is the initial ratio of carbon monoxide to carbon concentration. Secondly, the χ i are the ratios of the diffusivity in the porous layer to that in the interior and exterior. Thirdly, the mole fraction of carbon monoxide in the inert gas is denoted by C. Finally, M 2 and M 3 are the normalized mass differences of the silicon monoxide-carbon reaction, (2), and the Boudouard reaction, (7), respectively.

Dimensionless boundary conditions
Our far-field boundary conditions, (23), become On the interface r = 1 between the porous media and the exterior environment, (24) and (25), become Similarly, on the interface r = s 0 between the porous media and the gas cavity, we have On the surface of the quartz, r = s(t), boundary conditions (26)-(29), become where the constants measure how fast the carbon monoxide-quartz reaction, (6), proceeds relative to how quickly the gases diffuse on the surface of the quartz. Finally, the quartz consumption rate, (30), simply becomes is proportional to the ratio of the rate constants for the carbon monoxide-quartz reaction, (6), and the silicon monoxide-carbon reaction, (2). We note that not all of our dimensionless parameter groups are independent of each other. In particular, We consider the gas density, ρ 0 , the far-field carbon monoxide concentration, c 0 , and the density and molar mass of quartz to be constant. We thus consider the independent dimensionless parameters to be D i , s 0 , K, and S, and we determine Pe i and R i by (70).

Dimensionless initial conditions
The dimensionless versions of the initial conditions, from (35) and (36) are with

Quasi-steady regime
We estimate the values of the dimensionless parameters in Table 2. The wide range of some of the parameter values is primarily due to the uncertainty of the rate constants and radius of the pellets. Industrially, the most relevant limit is where the supplied concentration of carbon monoxide is low. Thus, in order to gain insight into the behaviour in an industrially relevant regime, we consider the regime, where δ C 1. Furthermore, we consider the distinguished limit where each of Pe i , s 0 , K, χ i , R i , and S are O(1); and D i is O δ −1 , to recover the richest possible dynamics. The resulting evolution of the gas concentrations are quasi-steady at leading order. The consumption of the carbon layer and quartz surface remain time dependent. There are three time-scales within our system. The short time-scale, the time-scale that we non-dimensionalized with, and the long time-scale. Equilibration of the gas concentrations happens on the short time-scale. Due to the density difference between the gases and the solids, the solids do not evolve on the short time-scale. On the scale we have chosen to nondimensionalize the system with, the carbon and quartz is consumed by the reactions and the gas dynamics are quasi-steady. On the long time-scale, all of the quartz will be consumed and thus, no silicon monoxide or carbon dioxide will be present in the system. On this scale, the system is static since no quartz remains and there is no silicon monoxide or carbon dioxide to react with any remaining carbon.
We take the limit in which C 1 and δ 1 in (41)-(44). In the exterior the leading-order quasi-steady model reads for r > 1; analogous equations hold in the gas cavity. The equations (48)-(50) for the evolution of carbon, silicon carbide, and the porosity remain the same in the limit δ 1. However, the conservation equations for the gases in the porous media Dimensionless diffusivity  (7) and (2) Dimensionless speed of (6) Relative volume change of (2) The typical values are computed using the 'value assumed' in Table 1 All the boundary conditions remain the same, except the boundary conditions (63)-(66), which simplify to

Analysis of the quasi-steady model
We analyse our quasi-steady model assuming D SiO = D CO = D CO 2 and χ SiO = χ CO = χ CO 2 for simplicity, and so Pe SiO = Pe CO = Pe CO 2 and R SiO = R CO = R CO 2 as a consequence. The analysis can be done in the general case, however, the manipulation is algebraically tedious and we do not expect the diffusivities to vary greatly among the three gas species. We drop the subscripts in the subsequent analysis. In the gas cavity and exterior of the pellet, we will solve explicitly the simple form of the partial differential equations governing the gas concentrations and, using these analytical expressions, we will reduce the problem to one posed in only the porous layer, with appropriate boundary conditions at r = s 0 and r = 1. A numerical solution is then required only for s 0 ≤ r ≤ 1. We first integrate the gas flow equations, (77) and (81), to find where we have used the boundary conditions (60) and (62). The velocities are thus determined by α(t), the function of integration interpreted as the time-dependent advective flux induced by the reactions. The α(t) is determined as a part of the solution, found later by (98). We substitute (86) and (88) into (74)-(76) to find that the concentration of gases in the gas cavity are for s(t) < r < s 0 , while in the exterior they are for r > 1, where the β i (t) and γ i (t) are to be determined. Eliminating β i (t) and γ i (t) using the boundary conditions (82)-(85) and (68), we find that, on the inner surface of the porous media, r = s 0 , the concentrations of the gases satisfy the conditions as well as the additional requirement that to determine α(t), which may be shown to be unique and non-negative by inverting (98). Similarly, on the exterior boundary, r = 1, we have We have reduced the model to a problem posed in just the porous layer, but we cannot reduce any further the system (48), (50), (68), and (78)-(80), with boundary conditions (95)-(101), and initial conditions (71) and (73). Hence, we seek a numerical solution using the method of lines using second-order central differences in space. This discretization turns the differential equations (78)-(80), and the boundary conditions (95)-(97) and (99)-(101) into a set of algebraic constraints, which combined with (48), (50), and (68), leaves us with a set of differential-algebraic equations (DAEs) to solve. We specify the resulting DAE system in Julia [32,33], and solve the system with the DAE solver IDA of SUNDIALS [34] using adaptive-step time integrations. Furthermore, we use the ModelingToolkit package [35] to efficiently and quickly compute the sparse Jacobian. We proceed to test the convergence in the usual way by refining the spacial grid and maximum tolerance.
We first consider the time evolution with the parameters given by the representative values in Table 2, shown in Fig. 4. We find that the carbon monoxide diffuses through the porous layer and reacts with the quartz core. As time goes on, the concentration  Table 2 rises since the reaction with the quartz slows as the carbon monoxide is consumed, and also because carbon monoxide is produced in the porous layer by reactions (2) and (7). The dynamics of the silicon monoxide is the converse of the carbon monoxide. The concentration of silicon monoxide is largest at t = 0, and decreases with time. Less silicon monoxide is produced on the surface of the quartz with time, and the carbon consumes the silicon monoxide by reaction (2). Since K = 1, the dynamics of carbon dioxide is identical to silicon monoxide. Overall, there is not a large variation in the concentration of silicon monoxide throughout the porous layer. As such, the carbon is consumed, and silicon carbide is produced, fairly uniformly. Since the concentration of silicon monoxide is largest at r = 0.8, the carbon is consumed, and the silicon carbide is produced, the fastest here. The radius of the quartz core decreases approximately linearly, since there is not a large variation in time of the concentration of carbon monoxide on r = s(t). However, since the concentration of carbon monoxide increases with time, the quartz core shrinks faster as time increases.
In order to better quantify the utilization of the carbon and the production of silicon carbide, we introduce two macro-scale metrics. The first metric is the carbon utilization fraction which is the normalized amount of carbon used, and is a measure of how much of the initial carbon reacts by either the silicon carbide reaction or the Boudouard reaction. The second metric is the silicon yield which is the amount of silicon carbide produced normalized by the initial quartz-it is the fraction of the initial quartz that has been converted to silicon carbide. In the example shown in Fig. 4, U C = 0.517 and Y SiC = 0.494 at the quartz consumption time. Only half of the carbon is used and the rest is wasted, and only half of the quartz is converted to silicon carbide and the rest is lost to silicon monoxide exiting the pellet. The relative strength of diffusion to the reactions within the porous layer is key to predicting the dynamics. As we saw, a significant amount of silicon may be lost to the transport of silicon monoxide out of the pellet. We now examine two limits of the simulations. We begin with the reaction-limited case, which we show in Fig. 5. In this case, we use the parameter values in Table 2, but with D = 10 5 . As we shall see later on, the choice D = 10 5 is sufficient to highlight the reaction-limited behaviour, despite Table 2 suggesting we choose D = 10 7 . Generally, we find the same qualitative behaviour as in the previous example. However, the spacial variation in the concentrations is reduced. Since diffusion dominates, the concentrations of carbon monoxide and silicon monoxide are uniform throughout space and time. Any silicon monoxide (or carbon dioxide) produced by the quartz reaction quickly diffuses through the porous layer and is lost, leading to slow consumption of carbon and production of silicon carbide. The losses are further reflected in the carbon utilization and silicon yield. At the quartz consumption time U C = 0.120 and Y SiC = 0.114.
We now examine the diffusion-limited case by using the typical values found in Table 2, but with D = 10 2 , with results shown in Fig. 6. Again, we choose D to show the diffusion-limited structure, not necessarily the minimum possible value. The concentration of carbon monoxide at the quartz interface is small, since the reaction dominates diffusion-almost all carbon monoxide that reaches the quartz surface reacts immediately. In turn, a considerable amount of silicon monoxide is produced. The concentration of silicon monoxide sharply decreases in the porous layer, since most of the silicon monoxide is consumed before it is able to diffuse to the exterior of the pellet. We find a rapid consumption of carbon and production of silicon carbide. In this case, we find a moving front-behind which, all the carbon has reacted-instead  Table 2, but with D = 10 5 of the relatively uniform consumption we saw in the previous examples. Owing to the moving front nature of the carbon, once all the carbon at a particular radius has been consumed, no more silicon carbide can be formed there. This can be seen with the highlighted curves in Figs. 6c and d. All the carbon has been consumed for r 0.9 and corresponds to where the silicon carbide curve begins to deviate from the maximum. The silicon carbide profile remains the same for r 0.9 for all later times because the carbon has been fully exhausted in this region, and no additional silicon carbide can be produced. In this example, the carbon utilization and silicon yield at the quartz consumption time are U C = 1.000 and Y SiC = 0.955. Since the reactions dominate diffusion, all the carbon is reacted, and almost all the initial quartz is converted into silicon carbide.  Table 2, but with D = 10 2 We now explore the parameter space, and use the carbon utilization, (102), and silicon yield, (103), at the quartz consumption time to characterize the effectiveness of the process.
In Fig. 7 we show the carbon utilization and silicon yield at the quartz consumption time for three characteristic planes within the parameter space. The base parameter values are D = 10 4 , s 0 = 0.8, K = 1, and S = 0.1, as in Table 2. We begin with the D-s 0 plane (Fig. 7a). For D < 10 3 we find a transition around s 0 ≈ 0.8 from maximal silicon yield to full carbon utilization. Since K = 1, the silicon carbide reaction, (2), and the Boudouard reaction, (7), dominate diffusion equally. Each mole of quartz is able to consume three moles of carbon-two moles by the silicon carbide reaction, and one mole by the Boudouard reaction. The initial carbon:quartz mole ratio is given by  Table 2, D = 10 4 , s 0 = 0.8, K = 1, and S = 0.1. We use slightly different parameter ranges than in Table 2 to better highlight the variations within the parameter space Thus, an initial mole ratio of 3 carbon:1 quartz corresponds to an initial quartz radius of s 0 = 0.793. Therefore, the pellet is quartz-limited for s 0 < 0.793, and carbon-limited for s 0 > 0.793, in the case of K = 1. We highlight the dependence of this critical s 0 value on the initial porosity in Fig. 8. Above the critical line the carbon:quartz ratio is less than 3, and the pellet is carbon-limited. Conversely, below the line the carbon:quartz ratio is greater than 3, and the pellet is quartz-limited. Hence, the silicon yield in Fig. 6f is less than 1 since the system is carbon-limited with s 0 = 0.8. We also find that, for large D values, the carbon utilization and silicon yield are reduced. The reduction is caused by diffusion dominating the reactions, silicon monoxide and carbon dioxide more easily escape the porous layer-slowing the carbon reactions, and reducing both the carbon utilization and silicon yield fractions. Moving to Fig. 7b, the D-K plane, we see a different behaviour. Generally, as K increases, the carbon utilization increases. By increasing K, the Boudouard reaction begins to dominate the silicon carbide reaction, and so more of the carbon is consumed. However, increasing K reduces the silicon yield, albeit weakly. As the rate of the Boudouard reaction increases, there is less available carbon for the silicon carbide reaction, and so the production of silicon carbide is slowed slightly. We see that increasing D decreases the carbon utilization and silicon yield. The carbon utilization behaviour is more interesting than in Fig. 7a, there is not a simple transition at a particular D value, instead the transition depends on K. Increasing K increases the speed of the Boudouard reaction, and so D must correspondingly increase in order for diffusion to dominate.
The final plane we will look at in detail is the S-K plane, as shown in Fig. 7c. We recall that S compares the speed of the quartz reaction to the silicon carbide reaction. The plane is split into four distinct regions depending on which reaction(s) dominate(s). In the lower left region, S 1 and K 1, and so the silicon carbide reaction is the fastest of the three reactions. Here, only about a third of the carbon is consumed and only about half of the quartz is converted to silicon carbide. Neither of the metrics approach one owing to the reasonably large D value-a non-negligible amount of the silicon monoxide diffuses through the carbon having not reacted. In the lower right portion of Fig. 7c, the reactions are dominated by the quartz reaction. The carbon utilization and silicon yield increase here owing to the higher concentration of silicon monoxide caused by the faster quartz reaction. The upper left region is instead dominated by the Boudouard reaction, (7). Any carbon monoxide created by the quartz reaction quickly reacts with the carbon via the Boudouard reaction, thus, the carbon utilization fraction is highest in this region. However, the faster consumption of the carbon prohibits the formation of silicon carbide, and the yield is correspondingly reduced. Finally, in the upper right region, the silicon carbide reaction is dominated by both the quartz reaction and the Boudouard reaction. The silicon yield is lowest in this region, since the silicon carbide reaction is the slowest of the three reactions. Interestingly, the carbon utilization fraction is not largest here. The dynamics within this region are diffusion-limited, and so we recover a moving front within the carbon layer, and the carbon layer shrinks from the inside. In turn, the capacity of the carbon to react with the silicon monoxide and the carbon dioxide is reduced. The gases can diffuse out of the pellet more easily-reducing the amount of carbon that is utilized.

Conclusions
In this paper, we have formulated and analysed a model for the evolution of a single quartz-carbon pellet under experimental exploration for use in a silicon furnace. The aim is to identify how the quartz can be efficiently converted to silicon carbide for subsequent processing to silicon. We have developed a spherically symmetric model describing the relevant chemical reactions and mass transfer within a single pellet. We then reduced the model to a quasi-steady regime, and solved the reduced model numerically. The quasi-steady model has four independent non-dimensional parameters, which allows us to explore the dynamics and parameter space relatively efficiently.
We found two distinct regions in the parameter space-a reaction-limited regime and a diffusion-limited regime. In the reaction-limited regime, most of the silicon monoxide produced on the surface of the quartz diffuses out of the pellet and is lost. However, in the diffusion-limited regime, most of the silicon monoxide is converted to silicon carbide owing to the reactions dominating the transport of the gas. We characterized the dynamics using two metrics, the carbon utilization fraction and the silicon yield at the time that the quartz is fully consumed. The carbon utilization fraction told us how much of the carbon had reacted, and corresponded to either the reaction-or diffusion-limited regime. The silicon yield told us how much of the initial quartz was converted to silicon carbide, and thus, how much silicon monoxide escaped from the pellet. We numerically found a critical ratio of initial quartz radius to pellet radius of about 0.8. Here, the pellet transitions from being carbon-limited to being quartzlimited. The 0.8 value is in agreement with a calculation based on the initial mole ratio. A diffusion-limited system maximizes the carbon utilization and silicon yield, since the silicon monoxide reacts with the carbon before diffusing out of the pellet.

Further work
We considered the small δ and C, large diffusion limit as each of these is physically relevant to silicon production. However, one possible continuation would be to investigate further sub-limits. Other distinguished limits are also with the ranges in Table 2 and may be instructive, for example C = O(δ) and D i = O(1). More experimental evidence could help reduce the wide parameter ranges and suggest particular sublimits to consider. In particular, reducing the ranges of the rate constants k 1 , k 2 , and k 3 . Asymptotic analysis of the sub-limits could provide insight into the mechanism for silicon carbide production.
Another consideration is the effect of temperature. Our model assumed the pellet is isothermal. However, this assumption is not likely to be the case in an industrialscale furnace. Taking temperature into account, the reaction rates would increase with temperature, probably following an Arrhenius equation. Incorporating the temperature dependence would allow us to find the optimal operating temperature. However, the relative speed of the reactions would be determined by the Arrhenius equation, and would not necessarily be independent. An instructive extension would be to find the optimal temperature to maximize the efficiency of the process.
Another extension is relaxing our assumption of neglecting gravity. We neglected gravity in order to maintain spherical symmetry, thus, our model and solution only depended on the radial distance. However, as the quartz core shrinks we would expect the quartz to remain in contact with the bottom of the carbon layer. Considering gravity, we would have an axi-symmetric model that is more representative of the underlying physics and geometry.
We directly wrote bulk reaction terms within the carbon layer, and estimated the effective rate constants. However, another approach would be to use the method of homogenization to derive the form and effective rate constant from the microstructure of the pores and surface reactions. Homogenization would allow us to encapsulate the pore-scale structure into our model on the scale of the porous layer. Furthermore, our single pellet model is the first stage in developing a model for the operation of a furnace containing many pellets. However, it would be computationally expensive to track the gas behaviour in every pellet. One alternative is to use the method of homogenization to derive an effective model that holds over the scale of the furnace and that captures all the microscale behaviour in the parameters. Homogenization has been used to describe similar complicated multiscale processes, for example, a grain model with applications to silicon processing [11,36], reactive decontamination [37], and reactive filtration [38].
Our ultimate goal is to use the insight gained from our mathematical model to assist with the design of new furnace protocols and ultimately reduce the amount of greenhouse gases emitted into the atmosphere.