Mathematical Modeling of the Production of Elastomers by Emulsion Polymerization in Trains of Continuous Reactors

A mechanistic model is proposed to describe the emulsion polymerization processes for the production of styrene–butadiene rubber (SBR) and acrylonitrile–butadiene rubber (NBR) elastomers in trains of continuous stirred tank reactors (CSTRs). A single model was used to describe both processes by choosing the proper physicochemical parameters of each system. Most of these parameters were taken from literature sources or estimated a priori; only one parameter (the entry rate coefficient) was used as an adjustable value to reproduce the kinetics (mainly conversion), and another parameter (the transfer to polymer rate coefficient) was used to fit the molecular weight distribution (MWD) experimental values from plant data. A 0-1-2 model for the number of particles and for the moments of the MWD was used to represent with more fidelity the compartmentalization effects. The model was based on approaches used in previous emulsion polymerization models published in the literature, with the premise of reaching a compromise between the level of detail, complexity, and practical value. The model outputs along the reactor train included conversion, remaining monomer composition, instantaneous and accumulated copolymer composition, the number of latex particles and particle diameter, polymerization rate, the average number of radicals per particle, average molecular weights, and the number of branches per chain.


Introduction
Two of the most important rubber products (copolymers) from the economical point of view-styrene-butadiene rubber (SBR) and acrylonitrile-butadiene rubber or nitrile rubber (NBR)-are produced via emulsion polymerization, mainly via the continuous cold process (5-15 • C) in trains of continuous stirred tank reactors (CSTRs). Some types of styrene-butadiene copolymer rubbers can also be produced via living anionic polymerization in solution, with a better control of the composition microstructure along the polymer chain, but the subject of this paper is the random copolymer produced via radical polymerization in emulsion. The world production of SBR was estimated to be around 8.3 × 10 6 ton/year in 2020 (based on [1]), and that of NBR was estimated to be 1.56 × 10 6 ton/year by 2023 [2], together representing roughly 3% of the total world polymer production (assuming a total world production of 300 × 10 6 ton/year).
The main use of SBR is in tires and automotive applications, while NBR is used in seals, gaskets, hoses, etc., where chemical resistance, higher strength, and lower gas permeability are needed. the use of pseudo-kinetic rate constants for the reaction kinetics, and the calculation of n by using the extended fraction expression [10]. It also used an average particle size and consisted of differential equations (dynamic model). As in the case of other models for continuous reactors, it was used to study start-up procedures and the control of polymer properties (copolymer composition, molecular weight, and particle size).
More recently, the Washington et al. model was used to study, in more detail, the effect of the monomer partitioning model, the radical desorption, and the monomer soluble impurities of several polymer properties in the batch and the CSTR train case. [30] The model was also applied to study dynamic aspects of a train of eight CSTRs and the effect of different operation policies on the off-spec material produced during transient operation. [31] Scott et al. also used this model to demonstrate the application of the Bayesian experimental design technique to maximize the information that can be experimentally obtained for an NBR process carried out in a train of CSTRs [32].
Given the similarities of the SBR and NBR processes, it was the goal of this work to develop and present a single model that is flexible enough to represent both processes with relatively small variations and to demonstrate that the same mathematical structure could be applied with minor adjustments to both processes. Some previous models have addressed only one of the two polymerization systems at a time (not necessarily for the continuous process but for batch or semibatch modes), although, in some cases, they used essentially the same model with adaptations. We believe that it is instructive to treat the two processes with the same model simultaneously in order to emphasize their similar underlying structures. Other goals of this work were to introduce improvements in some key aspects of the models (to be specified below) and to reach a good balance between the level of detail, the ease of implementation, and the predictive power of the model. To do this, a different model is proposed that is intended to study the effect of the configuration of the CSTR's train in the production (rate and quality) of both SBR and NBR elastomers. Configuration refers to the number of reactors in the train (in the range of 7-16), as well as the possible presence and position (reactor number) of side feed streams of one of the monomers or CTAs in the reactor train.
The goal of the study was to compare steady-state operations; nonetheless, a dynamic model was written in the form of ordinary differential equations (ODEs) since, in our experience, the numerical solution of this kind of system is more robust than that of the corresponding steady-state non-linear algebraic equation system. This selection was based on the experience of many years of our group working on the mathematical modeling of industrial emulsion polymerization processes.

Mathematical Model
In this section, the different components of the used mathematical model are presented and briefly discussed. Since some of these components were taken or adapted from existing models-only the main features of them are provided, and their details can be found in the given references and in the Supplementary Materials (SM).
The philosophy of the present work during the definition of the different components of the model was to use only the level of complexity deemed necessary to have a reliable description of each of the involved phenomena. In some cases, as in the calculation of n, which was based on a 0-1-2 particle model adapted from one by Vale and McKenna [33], this choice resulted in an increased level of detail compared with previous models. Additionally, the present model uses a more detailed calculation of the leading moments (0-2) of the MWD which, for consistency with the particle model, were also calculated using a 0-1-2 model based on a subset of the population balance model of Saldívar et al. [9] The model is not just a collection of existing pieces from other models, but the equations were re-derived and/or adapted as needed. In other cases, our choice was a simpler model than those in other works; this was the case, for example, of the monomer partitioning model, which is quite simple but still effective. A more detailed discussion of these differences is provided in the specific sections of the paper.
The used kinetic scheme is given in Table 1. Notice that the medium part of the table shows, in detail, the terminal model on which the pseudo-homopolymer approach was based (for propagation and chain transfer to monomers and CTAs). Diffusion-controlled termination using a single coefficient was used. The lower part of the table shows the kinetics in the pseudo-homopolymer form, which is particularly useful for the MWD (moments) calculations. Notice that the internal and terminal double bond reactions only occur on the butadiene units of the dead polymer, and this is reflected in the multiplication of the corresponding rate constant by F c 1,mo , which is the accumulated molar fraction of butadiene in the polymer. Table 1. Kinetics in the aqueous and particle phase.

Aqueous phase
Redox initiation Termination by combination and disproportionation (r,q = 1, . . . , ∞) 2 P r + P q ktc → D r+q P r + P q k td → D r + D q Chain transfer to polymer (q,r = 1,..,∞) 3 P q + D r rktrp → P r + D q Internal or terminal double bond polymerization (q,r = 1, . . . , ∞) 3 P q + D r F c 1,mo rk db → P r+q 1 A pseudo-homopolymer approach with apparent rate constants is used to represent the copolymerization kinetics based on the terminal model. 2 A single value k t = k tc + k td was used for the termination step, calculated as k t = √ k t11 k t22 while assuming diffusion-controlled termination, where k tjj , j = 1, 2, are the homotermination constants for monomers 1 and 2. 3 The transfer-to-polymer and double bond polymerization reactions are only considered for the MWD calculations.
The symbols used in Table 1 for the aqueous phase are: I w , initiator; Y r 1 and Y o 1 , reducing agents in the reduced and oxidized states, respectively; Y 2 , secondary reducing agent; R w , primary radicals in the aqueous phase; P w , polymeric radicals in the aqueous phase; and D w , the dead polymer in the aqueous phase. For the particle phase, the nomenclature is: P i , polymeric radicals of type i terminal unit; M j , type j monomer; D, the dead polymer; P r and D r , length (r) of live and dead polymers, respectively; M, total monomer (M = M 1 + M 2 ); and T, chain-transfer agent. Notice that P r = P r 1 + P r 2 , where P r i are r-length polymeric radicals of type i (i = 1,2).
The following equations are valid for any reactor in the train, except where noted.

Monomer Mass Balances
The remaining mass of monomer M m i (i = 1,2) Monomer bound in the polymer chain H m where M wi is the molecular weight of monomer i; w s, in is the inlet mass fraction of species s entering the reactor; F i,ma is the mass fraction of monomer i instantaneously formed in the reactor according to the Mayo-Lewis equation in mass form; F in and F out are the inlet and outlet total mass flows entering and exiting the reactor, respectively; M T is the total mass in the reactor; R p is the reaction rate, mol/(L s); and V w is the water volume in the reactor. Notice that the reaction rate is defined per liter of the aqueous phase. F ri,Mi are side feed streams (mass flow) of monomer i to reactor ri. The Mayo-Lewis equation in molar units is: where F i,mo (i = 1,2) is the mole fraction of i-monomer units instantaneously incorporated in the copolymer; f 1 (i = 1,2) is the mole fraction of monomer i remaining in the reaction site; and r i (i = 1,2) comprises the reactivity ratios of the copolymer system. The two forms (molar and mass) of the equation are related by a simple function f M that converts the molar fraction into the mass fraction: where M wi is the molecular weight of monomer i (i = 1,2). Notice that for practical reasons, mass, instead of molar units, is used for some of the balances, since the recipes are handled on a mass basis in industrial practice.

Rate of Polymerization
The polymerization rate can be calculated with the classical expression for emulsion polymerization: where [M p ] is the monomer concentration in particles, mol/L; n is the average number of radicals per particle; N A is the Avogadro number, mol −1 ; N p is the number of particles per L of water, L −1 ; and k p is the apparent rate constant for propagation, L (mol s) −1 , calculated as: where p i is the probability of type i radical in particle (i = 1,2) and ∅ i is the mole fraction of monomer i in particle (i = 1,2). The probabilities p i can be calculated as shown in [34].

Population Balance Equations for Particles and Calculation of R p and n
Defining F n,ri as the number of particles with n radicals (n = 0,1,2) per liter of water, L −1 , in the reactor r i in the train, population balance equations (PBEs) based on those postulated by Vale and McKenna [33] for a 0-1-2 system can be obtained. These equations are different from those of Vale and McKenna in two aspects: they do not consider the full-size distribution (partial derivatives with respect to size) but an average particle size instead that can increase with time and/or reactor. Another difference is that these equations include the inflow and outflow terms not present in the batch model of Vale and McKenna: dV w F 1,ri dt = ρF 0,ri − (ρ + k des )F 1,ri + (ρ + 2k des )F 2,ri + ρ mic M ic + Q in,ri F 1,ri−1 − Q out,ri F 1,ri (8) where ρ, k des , and k t are the radical entry, radical exit, and termination rate coefficients, respectively (see details and units below); v is the volume of an average particle, L; M ic is the micelle concentration L −1 ; Q in,ri and Q out,ri are the total in and out volumetric flows to and from reactor ri, respectively; . Notice that the radical entry and exit coefficients are particle size-dependent as detailed below. The micellar nucleation mechanism is included in Equation (8); in general, the radicals in the aqueous phase may enter micelles or particles in a competitive way.
A dimensionless version of these equations, used for the numerical solution of the model, can be found in the Supplementary Materials.
From the solution of Equations (7)-(9) it is possible to calculate n and N p from the following expressions:

Radical Entry and Exit Coefficients
The used specific expressions were taken from [9], but they were based on concepts previously proposed and generally accepted [35][36][37]. See the Supplementary Materials for specific details; Supplementary Material Equations (S1a)-(S1o).

Mass Balances of Species
In all the balance equations, w S p ,in represents the mass fraction of species S p in the reactor feed and M S p is the molecular weight of S p (except for M T that is total mass in the reactor). Differential equations representing the balances are written in molar units, except where noted. They are written for the initiator I w , reducing agents Y 1 and Y 2 , surfactant S, and chain transfer agent T r , and they are summarized below. More details on the assumptions and derivation of theses equations are included in the Supplementary Materials.
Initiator I w : Reducing agents Y 1 and Y 2 : For further details on the use of these equations, see Supplementary Material Equations (S1a-S1r) in the Supplementary Materials. Surfactant S: The surfactant balance is connected and solved simultaneously with the corresponding adsorption equilibrium (Langmuir isotherm) described by the Supplementary Material Equations (S1u)-(S1x) in the Supplementary Materials.
Given the total amount of surfactants by the mass balance, the partitioning of the surfactant adsorbed in particles (S a ), and free in water (S f ) can be defined; S f is given by a quadratic equation and includes the surfactant in solution and the surfactant in micelles. The micelle concentration is obtained by Equation (16) or (17): where [S] cmc is the critical micelle concentration of the surfactant (mol L −1 ), a em is the surface area per surfactant molecule, and r m is the radius of a micelle. Chain transfer agent (Tr): where F r,Tr is a side feed stream (mass flow) of CTAs to reactor ri. The partitioning equilibrium of the CTAs defines their concentration in the particles. This is based on the assumption that the component partitions between the particles and the aqueous phase according to the partition coefficient are defined as: where [Tr] p and [Tr] w are the concentrations of CTAs in the particle and aqueous phase, respectively. The simultaneous solution of the mass balance and the CTA partitioning is described in detail by Supplementary Material Equations (S1y)-(S1aa) in the Supplementary Materials.
Water (in mass units): Differential mass balance equations were also derived for primary radicals R w and polymeric radicals P w (aqueous phase), but the quasi-steady state approximation (QSSA) was assumed for these species, resulting in a simple quadratic equation for P w : where k mp is an entry coefficient for radicals in particles (m/s −1 ); r is the radius of an average particle; k mm is the entry coefficient for radicals in micelles; a m is the total micellar surface area; and M ic (L −1 ) is the micelle concentration.

Monomer Partitioning
The modeling approach used here was relatively simple but still effective. Madhuranthakam and Penlidis [30] concluded that the use of different models for monomer partitioning (e.g., partition coefficients and models based on equations of state) does not make a significant difference in the case of the more complex NBR system in which the water solubility of AN is important. For the simpler SBR case, an equipartition approach (similar proportions of the monomers in the droplets and particle phases) equivalent to that proposed by Dougherty [38], is used. Additionally, the approach used here allows for the a priori estimation of key parameters of the partitioning model based on the known physicochemical properties of the system components. The used modeling approach was based on the simple limit conversion for saturation, x sat , used by Gardon [39] as the conversion at which the monomer droplets disappear (end of interval 2). Up to this point, in the presence of excess monomers, the polymer particles are saturated with monomers and maintain a nearly constant concentration [M] p is calculated differently depending on the stage (interval) of the polymerization. The decision is made based on the equivalent conversion, x, defined as: where M 2w stands for the mass of monomer 2 in water. This latter quantity is nearly zero for SBR, but it is significant for NBR due to the partial solubility of AN in water. If where V p is the volume of this hypothetical particle, which can be calculated assuming volume additivity: where M wi = molecular weight of monomer i and ρ i and ρ pi (i = 1,2) are the densities of monomer i and homopolymer i, respectively. If x > x sat (mass conversion x, interval 3), [M] p is calculated assuming that all the remaining monomers, except for the possible presence of AN in the water phase (NBR case), are in the particles: where: and the superindex m indicates mass units. An approximate way of estimating the amount of AN in the water phase (M 2w ) in the NBR case is based on a partition coefficient, defined as the ratio of mass concentrations of monomer 2 in the particles and the water phase: As shown in Supplementary Materials Section 2 (see Supplementary Material Equation (S2d)), M 2w can be obtained by solving the following quadratic equation:

Total Mass Balance
Since there might be side-feed streams of monomers and/or CTAs along the train, a total mass balance per reactor (in mass units) is necessary: where M T,ri is the total mass present in reactor ri (reactors of different volumes can be considered along the train). Notice that in previous equations, the second subindex ri was omitted when it was clear that the balances were written for any reactor in the train. The same convention is used for the total mass flows entering and leaving the reactor (F in,ri and F out,ri , respectively). As an approximation, it is assumed that the time derivative is zero; that is, a constant mass is maintained in each reactor. In reality, a constant volume operation is used instead (level control); nonetheless, since the density variations are relatively small, this approximation introduces little error during the transient calculations and no error when the steady-state is achieved. In this way, the total output mass flow can be explicitly calculated from Equation (29) because the other terms are known (naturally, F in,ri = F out,ri−1 for ri > 1).

Molecular Weight Distribution
The quantities of two different distributions, one for the live polymer and the other for the dead polymer, are defined as follows: N r n = number of length-r radicals in particles having n radicals, per L of water. D r n = number of length-r inactive polymer chains radicals in particles having n radicals, per L of water.
For the assumptions and details of the derivation of the corresponding PBEs, the reader is referred to [9]. Though the derivations used here follow the general assumptions in [9], the model was adapted to comply with the 0-1-2 scheme. Some of the terms in the balances represent a shift in the classification of radicals in the presence of a given kinetic or physicochemical event.
A general balance for the live polymer r = 1, 2, 3 . . . and n = 1, 2, 3 . . . , in any reactor, where the identification of the reactor number (ri) is kept only for the inlet and outlet mass flows, is: where δ ( j = j 0 ) represents the Kronecker delta function, which is 1 when a given integer variable j = j 0 and is 0 elsewhere. For a dead polymer, the balance is, in principle, valid for r = 1, 2, 3 . . . and = 0, 1, 2, 3 . . .; the resulting PBE is: where w f w is the mass fraction of water at the inlet and P r , f n and D r , f n are the values of the distribution in the inlet stream for the live and dead polymers, respectively.
Since a 0-1-2 model was to be used, Equation (30) was applied only to n = 1,2, while Equation (31) was applied to n = 0,1,2. Then, the method of moments was applied to derive equations for the following moments: where µ K , n and λ K,n are the K-th moments of the live and dead polymers, respectively, in particles having n radicals. Only the first three moments (0, 1, and 2) of both populations are needed to calculate the number-and weight-average molecular weights. The final equations for these moments, six for the live polymer (Supplementary Material Equations (S3d)-(S3i)) and nine for the dead polymer (Supplementary Material Equations (S3j)-(S3r)), are provided in the Supplementary Materials, both in their original and dimensionless versions. The modeling approach followed here involved much more detail than previous modeling efforts since a full set of moments (0-2) had to be calculated for particles with 0, 1, and 2 radicals and then added together, taking compartmentalization effects into account more faithfully.

Number of Branches
The number of branches, B, per L of water in any reactor, can be estimated simply by the following ODE: In steady-state, this can be written as: or, equivalently:

Strategy of Solution
The model described in the previous section resulted in a first set of 13 ODEs per reactor (ri) associated to the states M 1 , M 2 , H 1 , H 2 , S, W, Tr, I w , Y 1 , Y 2 , F 0,ri , F 1,ri , and F 2,ri , as well as 15 additional ODEs per reactor for the moments of live (6) and dead polymers (9) (Supplementary Material Equations (S2d)-(S2r)). Note, however, that the equations for the live polymer moments µ 0,1 and µ 0,2 did not need to be solved, since equivalent information could be obtained from the solution of Equations (8) and (9) by realizing that µ 0, 1 ≡ F 1,ri and µ 0, 2 ≡ 2F 2,ri . Though the number of branches B could be dynamically calculated using Equation (34), it was decided to calculate its value only at the final steady-state by using Equation (36). In summary, the model resulted in a set of 26 ODEs per reactor. A summary of the ODEs and the most important auxiliary equations, which are the set of working equations implemented in the solution code, is included in Table 2.
To reduce the stiffness of the ODE system, the QSSA was applied to the remaining 4 live polymer moments. Simple explicit algebraic expressions were obtained for these variables, reducing the effective number of ODEs per reactor to 22.
Given the relative complexity of the equations (especially those related to the MWD) and to facilitate the programming and debugging steps of the implementation of the numerical solution of the ODEs, this task was implemented in two stages. In the first stage (computer program), only the first 13 ODEs (for each reactor) were integrated by using a given set of initial conditions in all the reactors in the chain until the steady-state was reached. The set of initial conditions was designed to reach the steady-state fast; linear profiles of the species along the reactor train, according to the approximate profiles expected at the steady-state, were used as the initial conditions. The evolution of the states was recorded at regular time intervals in a file, as illustrated in Figure 1. The variables in the first set of ODEs are independent of the MWD, and, therefore, the corresponding equations could be solved first. Once the solution of this set was obtained, the set of ODEs corresponding to the MWD was numerically solved in a dynamic way in a second computer program by using the dynamic information of the states obtained from the previous integration stage. To obtain more accuracy in the second stage, linear interpolation was used as needed to obtain values of the variables in the first set at integration times between the recorded time intervals (every~20 s). The numerical routine DDASSL (Differential-Algebraic System Solver) [40] was used in the two stages (programs) for the integration of the equations, which were coded in Fortran. DDASSL is capable of solving either a system of differential-algebraic equations or a system of pure ordinary differential equations, and, since the present system is numerically a system of ODEs (all the coupled algebraic equations are explicit), any ODE integrator for stiff equations could have been used. In this case, all the variables were declared as double-precision in Fortran, and the used numerical relative tolerance was 10 −7 -10 −6 . These numerical parameters allowed for safe and robust convergence in all cases. The execution of each stage took~1-2 s from time zero to the achievement of the steady-state in a standard laptop, even in the case of simulations for the maximum number of implemented reactors (16 in a train). The simulated time needed to achieve the steady-state was around 2-3 times the overall residence time of the reaction in the reactor train, which is usually close to 5 h. The implemented solution strategy turned out to be efficient, robust, and relatively easy to debug.

Parameter Values
The first problem to solve for the practical use of the developed computer programs was the selection of the physicochemical parameters of the simulated systems (SBR and NBR). Many of them were found in literature sources, but in several cases, the information was scarce or nonexistent. In the latter cases, most of the parameters were estimated a priori based on reasonable assumptions or indirect information; for example, the kinetic constants of the redox initiation system were estimated a priori based on data on the concentration of initiator at the last reactor exit. Notice though that these should not be considered adjustable parameters, since they were not fitted after comparison with the experimental data corresponding to the model outputs but estimated before the simulations were performed. A possible concern with this approach is how significant is the effect of these parameters on the final calculations and, therefore, if the error introduced with these estimations is significant. To answer this, the discussion can be better organized by classifying the prior-estimated parameters into three categories that are discussed next: -Kinetic parameters (k d1 , k d1 , f , C M 1 , C M 2 , C T 1 , C T 2 ): These were based on independent kinetic data (e.g., the kinetic constants of the redox initiation system, as discussed above) or assumed similar to values published for chemically similar systems; therefore, it is expected that little error could be introduced through these parameters. -Surfactant-related parameters (r m , a em ): These were also estimated by assuming similar values to those of chemically similar systems, which vary within relatively narrow ranges, so no significant errors are expected associated with these estimations. -Desorption-related coefficients (D wT , D pT , m d , m d T ): These parameters enter into the expressions of the desorption coefficient (see Supplementary Material Equations (S1h)-(S1i)) and do not have strong influence on the final responses, according to our experience. -Partition coefficients (K T , χ AN ): These parameters were inferred from known values of the solubilities of the respective components in water and organic media, so it is expected that the estimations used were not far from the real values. Indirect evidence of the adequacy of these estimations came from the behavior observed in the model for responses affected by these parameters, such as the copolymer composition (NBR case) and the molecular weights, which were close to experimental values.
Essentially, only one parameter was left as an adjustable value (entry rate coefficients to micelles and particles which were assumed equal, k mm = k mp ) to fit the kinetic data (mainly conversion), and a second parameter (transfer to polymer rate coefficient k trp ) was used to fit the MWD dispersity. This was justified since the value of the last parameter is usually model-dependent, as it is quite difficult to estimate it from independent experiments; see, for example, [28]. Since both the k trp and the k db affect the MWD dispersity and the number of branches B but it is practically impossible to discern their individual values based solely on the experimental values of the dispersity and B, it was decided to arbitrarily set k db = 0 and to fit only the value of k trp , which should then be considered as an effective value. In principle, it should be possible to measure by 13 C NMR techniques the number of trifunctional and tetrafunctional branches, mainly associated with transfer to polymer (one tertiary carbon) and internal double bond (two tertiary carbons) reactions, respectively, and from these values to independently obtain a rough estimate of both kinetic coefficients; however, this fell out of the scope of the present work. Regarding propagation and termination coefficients, it is well-known that the most reliable values for these parameters are those obtained by pulse laser polymerization (PLP) (when available), and there have been some reported for these coefficients for styrene and AN (see for example [41,42] for propagation coefficients); however, it was also interesting to compare the model results with kinetic parameters that had been used before by other groups working with emulsion systems, in particular for styrene polymerization, [43] or specifically for NBR polymerization, [24,28], even though they were not determined by PLP. Both types of parameters were tested in the present model; however, for comparison with previous work, the values shown in the simulations were not necessarily obtained with PLP-determined parameters (see details in Table 3). Nonetheless, the sensitivity of the model outputs to the used parameter type (determination method) was quite low. For example, for the SBR system, by changing the value of the propagation coefficient k p,22 for styrene from the non-PLP parameter value (151 L mol −1 s −1 at 10 • C) [43] to the PLP parameter value (43 L mol −1 s −1 at 10 • C) [41] for simulated conditions very similar to those in Table 4, there were only minor changes in the final conversion in a 10-reactor train (from 0.740 to 0.732,~1%) and in the number of particles (from 2.74 × 10 18 to 2.76 × 10 18 L −1 , less than 0.5%), as well as even smaller changes in all the other outputs. This low sensitivity can be explained in part due to the relatively low amount of styrene used in the SBR formulation (~17% on a molar basis), but other mechanistic features of this kind of system could also play a role. For the NBR system, a similar situation occurs. A change of the propagation coefficient value k p,22 for AN from the non-PLP parameter (6630 L mol −1 s −1 at 10 • C) [24] to the PLP parameter value (2604 L mol −1 s −1 at 10 • C) [42] for an initial monomer feed of f 1 = 0.68 (wt. basis (without additional side stream of AN)), results in changes smaller than 0.2% in all the model outputs. Unfortunately, for butadiene, the available PLP-determined expression for k p,11 [44] is not applicable in the temperature range of interest for these processes; the same is true for the termination coefficient for butadiene, k t,11 , but in this case, it was decided to extrapolate the Arrhenius expression in [45] out of its valid temperature range because the activation energy and, therefore, the temperature effect is much smaller for the termination than for the propagation step. Table 3 contains the values or range of values used in the simulations. For proprietary reasons, some of the details of the formulation and specific parameters used are not disclosed.  In all kinetic constants k x where an Arrhenius expression was available, the used notation is: with R = 1.987 cal mol −1 K −1 and the temperature T in K. 2 Values used in all the simulations in figures. 3 Only the used values that are different from those of the SBR system are listed.
Based on ferrous sulfate. 2 Based on sodium formaldehyde sulfoxylate. Figure 2 shows a comparison of the model predictions with plant data for conversion (18 experimental points) in the SBR case for a 10-reactor train. The steady-state operation of each reactor (except for reactor 2 (R2) in the train) was sampled at two different times (data labeled as experiments 1 and 2); slight temperature variations (± 1 C) were recorded for the temperatures of the reactors at these different times, except for reactor 10, in which the temperature difference was +3 C for experiment 2. Details of the used temperature profile, as well as the formulation data used in the plant, are provided in Table 4. A similar formulation has been published before for the NBR case [29,32]. The agreement between model and experiment was quite reasonable using a single value of k mm = k mp .  Table 4. Figure 3 shows the model predictions with a few experimental data from the plant available for the number average molecular weight M n , the MWD dispersity, and the accumulated copolymer composition F 1 . The copolymer composition was measured by refractive index (ASTM D5775-95(2019)), while M n and the MWD dispersity were measured by size exclusion chromatography (estimating absolute values from those relative to a polystyrene standard). The reactivity ratios were slightly adjusted (~10%) with respect to those reported in the open literature to get a better agreement with the copolymer composition data (Table 3), but this was justified since the literature data were obtained at a different temperature. The agreement of the model with the experimental data was reasonable but unfortunately limited to the few available experimental data points. Figures 4-6 show the effect of variation of surfactant on the most important outputs and process parameters with respect to the base case, which corresponded to the conditions of Experiment 1 and included a side stream of CTAs in R5. The predicted effects were all expected and easy to explain. As the surfactant increased with respect to the base value, more particles were produced ( Figure 5A) and the reaction rate was higher ( Figure 5C), as was the conversion (Figure 4A). At a higher conversion, the composition drift increased, as shown in Figure 4B-D. A higher number of particles resulted in smaller particle diameters ( Figure 5D), although the effect is barely perceptible in the plot. Finally, higher reaction rates resulted in higher molecular weights and dispersities. The jump observed in R p in R6 was due to the abrupt increase in temperature from R5 to R6 (+10 C). Figure 5B shows the value of the average number of radicals per particle ( n) according to the model. It must be highlighted that some of the previous models ( [3] and partially in [4,5]) assumed a constant value of n = 0.5 (Smith-Ewart case II kinetics), but, as seen in the figure, this may have deviated from the presumably true value by 20-30%, thus having a direct impact on the estimation of the polymerization rate (see Equation (5)).  Table 3, and the reaction conditions were similar to those in Table 4.  , reaction rate R p (C), and particle diameter D p (D) along the reactor train with respect to the base case in SBR production. Figure 6 shows how an increase in surfactant concentration also caused an increase in the average molecular weight (Figure 6A,B). This was due to an increased number of particles that delayed the entry of a second radical to the particle (more competition for radicals), extending its life and, therefore, its molecular weight. The increased competition for radicals also resulted in an increased dispersion in the probability of capturing a second radical, leading to increased MWD dispersity ( Figure 6C). Notice that this explanation has also been put forward by other researchers, albeit in the frame of miniemulsion polymerization. [52,53] Figure 7 shows the evolution of the number of particles with zero, one, and two radicals for the simulated SBR base case. At the first reactor, due to the relative abundance of initiator radicals, the number of particles with one radical (F 1 ) was higher than that with zero radicals (F 0 ); however, down the train, this situation was reversed as the initiator was depleted. From the plots, it is also clear that the number of particles with two radicals (F 2 ) was negligible in all the reactors, pointing to the conclusion that a 0-1 model should suffice to represent these systems. Notice that this behavior does not mean that an average value of n = 0.5 would provide a good representation of the system as assumed in early models; it rather means that n < 0.5. However, the whole situation could change if a gel effect was included in the model.  Figure 8 shows the effect of the variation of the amount of CTAs on conversion, average molecular weights, and MWD dispersity. The effect was almost negligible on the conversion since the desorption of CTAs was not significant; however, the effect on the molecular weight averages and dispersity was as expected, with large increases in these quantities as the amount of CTAs decreased. The plot for a 45% decrease of CTAs is not shown since the calculation diverged for the last reactor as the system reached the gelation point (the second moment of the MWD tended to infinity). Finally, Figure 9 shows the effect of the side stream of CTAs on the MWD. As mentioned above, in the base case, a (correction) side stream of CTAs was fed to R5. Figure 9 compares the behavior of the train with and without the intermediate CTA feed. As previously discussed, the effect was negligible for the progress of the conversion along the train ( Figure 9A), but the final M n and the dispersity increased from~70,000 to~90,000 and from~7 to~9, respectively, when the CTA side stream was suppressed. Similarly, the final number of branches increased from~1 to~1.2 in the absence of the side stream. Clearly, the intermediate CTA stream avoided an excessive increase in the molecular weight and the eventual formation of gel in the last reactors of the train that operate at high conversions. The effects of the changes in the redox initiator system were milder than those found for the surfactant and CTAs and are shown in the Supplementary Materials ( Figure S1). Similar behavior was found for the NBR system upon changes in the surfactant, CTA, and initiator. The effects of changes in the surfactant for an NBR system are shown in Figures 10-12 for the same important outputs and process parameters shown for the SBR system in Figures 4-6. The variations of the surfactant were made with respect to a base case, which corresponded to a composition feed of f 1 = 0.68 (wt. basis) and weight feed ratios of 1.25/0.082/0.35/100 for S/I w /CTAs/monomers. The base case also included a side stream of an additional 20% flow rate of CTAs at R6. The used temperature was 10 • C in all reactors except in R1 (12 • C) and R10 (15 • C). The effects were very similar to those found for the SBR system, with a few exceptions. In Figure 10, the butadiene composition of the remaining monomer f 1 , as well as those of the copolymer (accumulated and instantaneous (F 1 and F 1,inst , respectively)), increased instead of decreasing like in the SBR case. This was expected given the differences in the reactivity ratios between the two systems. In this region of the monomer feed composition (above the azeotropic point for NBR) in the NBR case, acrylonitrile was consumed faster than butadiene (see the Mayo-Lewis plot for this system in Figure S2), and the composition of the remaining monomer therefore increased in butadiene with conversion, also leading to a richer butadiene composition in the copolymer. Figure 10. Effect of variation of surfactant (from −45% to +45%) on conversion (A), remaining monomer feed composition f 1 (B), instantaneous copolymer composition (F 1,inst ) (C), and accumulated copolymer composition F 1 (D) along the reactor train with respect to the base case in NBR production. All compositions are weight fractions. Figure 11. Effect of variation of surfactant (from −45% to +45%) on number of particles N p (A), average number of radicals n (B), reaction rate R p (C), and particle diameter D p (D) along the reactor train with respect to the base case in NBR production. The effects of the surfactant changes on N p , R p , n, and D p shown in Figure 11 for the NBR system were qualitatively very similar to those observed and already discussed for the SBR case ( Figure 5). The only minor difference was that the effects were smaller for n in the NBR than in the SBR case. Another difference is that in the NBR case, there was not a jump in the reaction rate as in the SBR case; however, this was simply because in the former case, there was not a jump in an intermediate reactor temperature as in the latter. Figure 12 shows similar effects of the surfactant concentration for the NBR system compared to those in the SBR system ( Figure 6). A difference was that in the NBR case, the molecular weights reached a peak in R5 to decrease in R6 due to the side stream of CTAs. This effect was less important in the SBR base case, but in that system, the higher temperatures in the last reactors also significantly contributed to higher molecular weights. Additionally, the molecular weights and dispersities were lower for NBR compared to those exhibited by the SBR system.

Model Validation
The effects of changes in the CTA feed with respect to the NBR base case are shown in Figure 13, which is comparable to Figure 8 for the SBR system. Again, the behavior of both systems was qualitatively similar with the only exception of the peak exhibited by the molecular weights in R5 for the NBR case due to the side stream of CTAs in R6, as in Figure 12, and the different temperature profiles used for the SBR and the NBR cases. More interesting behavior for the NBR system was related to the drift in copolymer composition along the train. Given the reactivity ratios of this system (0.30 and 0.04), the system exhibited an azeotropic composition at f 1,az = F 1,az = 0.578 (molar basis, 0.582 weight basis) (see the F 1 vs. f 1 Mayo-Lewis curve in the Supplementary Materials, Figure S2) Therefore, composition drift depends on what side of the azeotrope the composition of the feed is located. If the feed composition is above f 1,az , as in most commercial grades of NBR, without correction, the instantaneous product composition will initially deviate towards copolymers with a lower Bd composition than the feed (R1); however, as the reaction mixture proceeds along the reactor train, the Bd-enriched environment will tend to increase the Bd content in the copolymer. This can be seen in Figure 14 for a product grade in which f 1 in the feed was 0.68 (wt. basis). Since composition drift may lead to a product with heterogeneous composition or a product deviated from the target average composition, it is common to introduce a side feed stream of AN in an intermediate reactor of the train to "correct" or attenuate the composition drift. Figure 10 also shows the composition profile when an additional AN side stream (20% of the initial feed) was fed to reactor R7. Without the correction, the final F 1 was above 0.66; on the other hand, with the side feed stream at R7, the final F 1 was slightly below 0.65. A second example is shown in Figure 15; in this case, the composition of the feed was f 1 = 0.5 on a wt. basis, a composition which was located on the other side of the azeotrope (see the F 1 vs. f 1 curve in the Supplementary Materials). The behavior of this system was the opposite of that in the previous example, exhibiting a composition drift that initially produced copolymers with a higher Bd composition than the feed (R1, F 1~0 .58) and a decreasing Bd content in the copolymer along the reactor train to end up with an F 1 slightly above 0.56, although the instantaneous F 1 fell below 0.54 at R10. Figure 15 also shows how this drift could be attenuated by feeding a side stream of Bd (20% of the initial feed). This correction worked very well, resulting in a final average composition near the target (F 1~0 .58) while avoiding large deviations in the instantaneous composition produced in each reactor (better composition homogeneity). Figure 14. Effect of the operation of the reactor train for NBR production with (corrected) and without a side stream of additional AN in R7 (20% with respect to the AN feed in R1) on conversion (A) and instantaneous butadiene copolymer composition (F 1,inst ) and accumulated butadiene copolymer composition (F 1 ) (B). The initial monomer feed was f 1 = 0.68 (wt. basis). Figure 15. Effect of the operation of the reactor train for NBR production with (corrected) and without a side stream of additional Bd in R7 (20% with respect to the Bd feed in R1) on conversion (A) and instantaneous butadiene copolymer composition (F 1,inst ) and accumulated butadiene copolymer composition (F 1 ) (B). The initial monomer feed is f 1 = 0.5 (wt. basis).

Conclusions
A single mechanistic model was built and used to describe the production of both SBR and NBR elastomers in trains of emulsion polymerization CSTRs. The aim was to balance the level of complexity of the model with its capability of describing the effects observed at the plant on the process parameters in a quantitative way and the product quality upon variable changes. Some novel aspects of this model compared to previous models are summarized below, along with additional comments: - The radical compartmentalization was considered in more detail, both in the particle kinetics description and in the calculation of the moments of the MWD, and for the first time, a 0-1-2 model was used to describe these two important aspects in these processes. It was found that at least a 0-1 model should be used for an accurate description of the compartmentalization effects, especially for the first reactors in the train. A similar conclusion was recently put forward by Marien et al. [52] using stochastic simulation, although for a general miniemulsion system. - The monomer partitioning model was as simple as possible, and its parameters could be estimated a priori based on published physicochemical data. As mentioned by Madhuranthakam and Penlidis,[30] the use of more sophisticated thermodynamic models for monomer partitioning seems to be unnecessary for this kind of system. -A single modeling framework was presented in a unified form for the SBR and NBR systems. Though previous models have been applied with adaptations to both systems in some series of papers, it is believed that by presenting them together, it is easier to appreciate the similarities and differences in both systems.
The model is presently used in the plant to design changes of configurations of the reactor train, such as the number of reactors in the train (up to 16 reactors in the existing facility), the reactor volume, the location and flows of side streams of CTAs and monomers, and grade recipes, among other variables, to maximize the plant operation versatility maintaining high productivity and product quality. The model is quite flexible, and it has been implemented with the help of a friendly graphical interphase. Current work involves the generation of product quality correlations that link model outputs, e.g., M n and M w , and number of branches, with product performance parameters, such as Mooney viscosity.