NC State
Kerekes, R., McDonald, D., and Zhao, J. (2020). "Perspectives on deriving mathematical models in pulp and paper science," BioRes. 15(4), 7319-7329.


Mathematical modelling is a powerful tool in science. Causal mathematical models based on a clear picture of how key variables interact enable a deeper understanding of a given situation and provide reliable predictions. This is a classic approach in science. Unfortunately, this approach is declining in pulp and paper-related research in favour of simply reporting experimental data. The lack of a framework provided by a model diminishes the value of much experimental work. Therefore, the increased use of mathematical models is encouraged, and this approach is illustrated via several practical examples from our work.

Download PDF

Full Article

Perspectives on Deriving Mathematical Models in Pulp and Paper Science

Richard Kerekes,a David McDonald,b,* and Joe Zhao c,d

Mathematical modelling is a powerful tool in science. Causal mathematical models based on a clear picture of how key variables interact enable a deeper understanding of a given situation and provide reliable predictions. This is a classic approach in science. Unfortunately, this approach is declining in pulp and paper-related research in favour of simply reporting experimental data. The lack of a framework provided by a model diminishes the value of much experimental work. Therefore, the increased use of mathematical models is encouraged, and this approach is illustrated via several practical examples from our work.

Keywords: Models; Mathematical models; Papermaking; Research methodology

Contact information: a: Pulp and Paper Centre, 2385 E Mall, University of British Columbia, Vancouver, BC, Canada V6T 1Z4; b: JDMcD Consulting Inc., 97 Rue Kerr, Vaudreuil-Dorion, Quebec, Canada J7V 0G1; c: Tri-Y Environmental Research Institute, 2655 Lillooet St. Vancouver, BC, Canada V5M 4P7; d: State Key Laboratory of Biobased Materials and Papermaking, Qilu University of Technology, Changqing, Jinan, Shandong, P. R. China, 250353;

* Corresponding author:


Although there are differing ways of carrying out scientific research, a classic approach is one of visualizing the underlying mechanism of a problem and then describing it via a mathematical model. The resulting equations can provide quantitative predictions of complicated, non-linear processes over a wide range of operating conditions and parameters. As Wigner (1960) observed, it is amazing that mathematics is so effective in describing nature. From a practical standpoint, the models permit optimization of existing processes, solutions for performance issues, or evaluation of new technologies. Despite this value, this approach is declining in pulp and paper research, thereby diminishing the value of much experimental work. To encourage researchers to employ more mathematical modelling, several examples from our research are presented to illustrate how some past models were developed.

Pulp Fibre Flocculation

Pulp fibres aggregate into mass concentrations called flocs that produce non-uniform mass distribution in paper and a yield strength that affects flow rheology in process equipment. Given these effects, there was a need to quantify flocculation in a simple way.

Early work by Mason (1950) showed that fibre flocculation was primarily mechanical in nature due to inter-fibre contact. He found the condition for incipient flocculation to be a critical concentration, defined as one fibre in a volume swept out by the length of one fibre. This visual image suggested a broader picture of multiple fibres in this volume to characterize flocculation in various states (Kerekes et al. 1985). This concept was further explored in the form of a Crowding Number, N. For pulp fibres, this number can be calculated from Eq. 1,


where C is suspension consistency (%), l is fibre length (m), and is fibre coarseness (kg/m) (Kerekes and Schell 1992). Mason’s critical concentration is = 1.

Subsequent work by Meyer and Wahren (1964) showed that fibre suspensions at higher concentrations developed a yield stress when fibres became locked in bent configurations in the network, and friction restrained movement between them. This condition required fibre contact with three other fibres, which Soszynski and Kerekes (1988) showed to occur at approximately = 60.

In a later study of the gravity settling of fibres, Martinez et al. (2001) identified a third critical value, = 16. The physical importance of this relative to N = 1 and N = 60 was unclear. In later work, by making comparisons with the Effective Medium and Percolation Theories of fibre networks, Celzard et al. (2009) found that = 16 corresponded to a “connectivity threshold” whereas = 60 corresponded to a “rigidity threshold”. This finding implied that, at = 16, each fibre is in contact with two other fibres. These three regimes of behavior are illustrated in Fig 1.

Fig. 1. Illustration of fibre contacts for various values of N: (a) Occasional contact (= 1); (b) Connectivity threshold (N = 16); (c) Rigidity threshold (N = 60)

Thus, fibre suspension behaviour changes dramatically from occasional contact among fibres to a network possessing strength in the relatively small range of 1 < N < 60. As a result, the Crowding Number can be employed to set desirable consistencies for various processes. For example, commercial papermaking requires the highest fibre concentration that can be readily dispersed. For this reason, headbox consistencies are typically set in the range of 20 < N < 60. Thus, lower consistencies are needed for longer fibres. The same holds true for hydrocyclone pulp cleaners. It is noteworthy that standard handsheet making, which requires minimal or no flocculation, employs a consistency of approximately 0.02%, which corresponds to approximately = 1 for a 1.5-mm fibre.

In summary, an equation based on a simple concept enables the setting of appropriate consistencies for differing pulps in pulping and papermaking processes.


Separating exceptionally large particles from small ones by screening is a trivial problem. However, separating nearly identical particles with high aspect ratios is complex. This is the case of fine screening pulp, such as the separation of shives from fibres. For high throughput, aperture sizes must be larger than the minimum dimension of either accept or reject components. This requires probability screening, in which separation is governed by factors such as particle flexibility rather than physical obstruction of the reject component on any of its dimensions.

As both accept and reject particles may pass through apertures, the efficiency of probability screens (E) is highly dependent on the amount of feed flow forced through the apertures. The fraction of the feed flow that is not forced through the aperture is the reject ratio (R). Separation efficiency (E) is described by performance equations governed by and a screen parameter (α). The most common performance equation used in the industry is the Nelson (1981) equation shown below (Eq. 2):


An important need in screening was to predict the possibility of plugging due to thickening of the feed flow along the screen plate. When used to develop a thickening factor, the Nelson (1981) equation failed, but the reason for this was not clear. This suggested exploring a rarely used screen performance equation, the Kubat-Steenberg equation (Kubat and Steenberg 1955), which is shown below with screen parameter β (Eq 3):


Equations 2 and 3 have the same limits, i.e.= 0 at R = 0, and = 1 at R = 1In addition, Eqs. 2 and 3 have virtually the same curve shapes between these limits, and the screen parameters are linked in the form: α = 1 – β.

Fig. 2. The plug flow model correctly predicts high thickening at low reject rates, whereas the mixed-flow model does not (Gooding and Kerekes 1992)

To determine how Eqs. 2 and 3 are different and the importance of their differences, Gooding and Kerekes (1989) derived both equations from first principles. They found that the Nelson (1981) equation assumed “fully mixed flow” in the feed flow, whereas the Kubat-Steenberg equation assumed “plug flow” in this zone. The reason for the failure of the Nelson equation to predict plugging became obvious. If the feed flow is “fully mixed,” it cannot thicken along the screen plate. In contrast, “plug flow” permits thickening. A successful thickening factor was developed based on the plug-flow model, as shown in Fig. 2 (Gooding and Kerekes 1992).

A practical outcome of this work was a recommendation for mill simulation programs to employ the plug flow model (Kubat-Steenberg equation) in place of the Nelson equation to accurately account for the influence of the reject rate on reject thickening. In summary, successful modelling requires a clear picture of the physical problem and assumptions underlying the equations employed to describe it.


During the 1970’s and 1980’s, the use of twin-wire forming for high-speed papermaking grew rapidly. The most common type of twin-wire forming was a blade former, in which pulp trapped between two fabrics (screens) was drained when passing over stationary blades. Such formers resulted in excellent formation but poor retention of fine material compared to roll formers; the cause was unknown. Indeed, there was no clear understanding of how drainage occurred. This was the outstanding problem in papermaking circa 1990.

Fig. 3. Illustration of blade formers

Zhao and Kerekes (1995) examined the hydrodynamics of blade gap formers. The problem was visualized as a pressure build-up between two wires as they “pinched” when wrapping around a blade, as shown in Fig. 4.

Fig. 4. Model of twin-wires passing over an infinitely thin blade

Employing appropriate simplifications and approximations, the hydrodynamic analysis yielded a pressure distribution in the direction of motion, p(x), which is shown by Eq. 4,



where U0 is the fabric speed (m/s), H0 is the upstream gap size between fabrics (m), T is the fabric tension (N/m), ρ is the suspension density (kg/m3), k2 is the permeability of the formed portion of the web (m/s/Pa). Angles α1 and α2 (rad) are shown in Fig. 4.

Equation 4 showed that pressure (p(x)) increased sharply in a zone upstream of the blade edge. This pressure gradient caused a relative velocity between the pulp and wire in the machine direction, which was shown later to account for formation improvement. The accompanying large peak pressure accounted for low retention.

The model provided several important insights. For example, some blade formers produced an unexpected two-sidedness in paper structure. The finding that pressure formed upstream of edges meant that, if the blade was wide, water drained through both wires at the leading edge but only on one side of the trailing edge. Thus, only half the drainage in this “twin-wire former” took place between two wires. Further work showed that reducing the blade width to approximately 10 mm caused the two pulses to merge into one pulse upstream on the leading edge, thereby resulting in two-sided drainage (Green et al. 1998). Subsequently, narrow blades were employed in these formers. In summary, a mathematical model gave new insights that led to better operations.


The removal of water by wet pressing is a key step in papermaking. It is a complicated non-linear process that depends on many variables. A mathematical model is needed to optimize the process. Many models have been proposed, but they were either not practical or lacked the accuracy to evaluate optimization strategies and new technologies (McDonald and Kerekes 2017a). The “picture” on which most models were based was one of parallel pressures: one pressure on fibres and another pressure on water. As the web was dewatered, pressure shifted from water to fibres, which provided a structural support that shielded the remaining water from pressure. Thus, flow decreased and stopped because the pressure on water decreased and stopped.

In the cited study, an alternative picture was proposed, in which one pressure acts on both fibres and water. Decreasing permeability and not decreasing pressure causes flow to decrease during dewatering. This was justified on the basis that wet pulp fibres lack the structural integrity to shield water from pressure in a web that has more water than fibre, even at the end of pressing. Pressure is transmitted through fibres to the water within them. Consequently, the rate of water removal decreases with time because water is first expelled from large pores, then small pores, then eventually ceases altogether when water has no pathway to the fibre surface or is held there by surface tension. This model was entitled the Decreasing Permeability Model (DPM) (Kerekes and McDonald 1991; McDonald and Kerekes 1991). This picture was largely rejected at the time of its publication.

The above picture was modeled by Darcy’s law having a permeability dependent on a power function of moisture content (Kerekes and McDonald 1991; McDonald et al. 2000). For “flow-controlled” pressing, this gave Eq. 5,


where I is the press impulse (nip load/speed) (kPa.s), W is the basis weight (kg/m2), ν is the kinematic viscosity (m2/s), and the furnish is characterized by the specific permeability (A) (g/m) and the compressibility factor (n).

Equation 5 gave excellent fits to data from commercial machines and pilot presses (McDonald and Kerekes 1991; McDonald et al. 2000). In addition, it predicted the rule-of-thumb that a 10 °C change in temperature increases the pressed solids by 1% through the kinematic viscosity (ν). Over time, the DPM was used to evaluate optimization strategies and press design using the appropriate furnish dependent coefficients (A and n) for different grades.

As experience with the DPM grew, some applications that involved factors beyond those in Eq. 5 were encountered. One such case was lightweight papers, in which rewetting contributed substantially to the final moisture, and pressing took place to near equilibrium conditions. The latter made dewatering “pressure controlled” rather than “flow controlled,” as modelled in Eq. 5. To account for these factors, the DPM was extended to include a rewet term (R) (McDonald and Kerekes 1995) and an equilibrium term (me) (Kerekes et al. 2013). The full equation is shown in Eq. 6, which covers the full range of wet pressing, as shown in Fig. 5.


Among other applications, the DPM has been employed to investigate the factors that control the limit of pressing. Based on Eq. 6, McDonald and Kerekes (2017b) showed that rewet and equilibrium moisture controlled pressing to high dryness. Further, they showed that felt design has a major effect on rewet (McDonald et al. 2013) and that rewet occurs in two stages: flow from the felt back into paper and film splitting when these separate (McDonald and Kerekes 2018).

The factors that determine equilibrium moisture content me were still unknown. It was hypothesized that me is directly related to the surface tension of water in the fibre pores (Kerekes and McDonald 2020). This led to the development of an equation that relates the ultimate possible dryness in pressing to the structure of the fibre.

Fig. 5. The DPM describes the full range of wet pressing (McDonald and Kerekes 2017a)

In summary, a mathematical model built on a more realistic picture than one in common use at the time led to useful predictions and insights into wet pressing. The path to the final equation evolved over time as experience was gained for various applications.


In some machine calenders, paper may be separated from the roll surface by pockets of air between the paper and the roll (Fig. 6). This bubbling typically forms on the wrap preceding the bottom nip. It is most pronounced for impermeable papers that are being heavily calendered at high speed. The bubbles can cause paper to fold over, which leads to creases or cuts that can cause breaks during reeling, winding, or converting.

Fig. 6. Bubble position in the cross direction (Hamel et al. 1999)

The prevailing theory on the cause of bubbles was that cross-direction sheet expansion was the key factor. However, this could not explain why increased machine speed or decreased paper permeability increased the onset of bubbling. The problem defied a mechanistic description until it was proposed that air expelled from the ingoing nip was responsible for bubble formation (Hamel et al. 1999).

The mechanism of bubble creation was pictured as follows. Air in paper is expelled when paper is compressed in a calender nip. As the air cannot pass though the nip, it must flow back out but can only do so on the ingoing side. Air expelled from the side in contact with the roll must flow back through the paper thickness or accumulate between the paper and the roll. The accumulation creates a bubble whose size depends on the rate of air expulsion from the nip and the rate of flow through the paper. The latter depends on the permeability of the paper and the pressure built up in the bubble. This can be solved analytically by imaging the equivalent problem of a calender roll rolling against paper on an infinite plane as shown in Fig. 7.

Fig. 7. Simplified bubble geometry (Hamel et al. 1999)

Hamel et al. (1999) modelled this mechanism by employing a force balance that assumed the paper had negligible stiffness. This led to Eq. 7,


where Ymax is the maximum bubble height (cm), l is the bubble length (cm) in the machine direction, ϵn is the in-nip compressive strain, zi is the thickness (μm) of the paper before the nip, S is the machine speed (m/s), K is the air permeability measured under applied pressure (m3/N.s), and T is the sheet tension per unit width (N/m).

Equation 7 explains several observations. Calendering decreases paper permeability and web tension in successive nips due to inelastic expansion in the machine direction. Therefore, bubbling can be reduced by several means, such as increasing the draw (tension) between the last dryer and the calender, pre-calendering, increasing calendering in the upper nips, increasing paper permeability, or improving cross-machine profiles (basis weight and moisture affect tension). The effect of these factors on bubbling is quantified in Eq. 7.


Crepe wrinkles are narrow machine-direction folds produced within paper rolls during reeling and winding (Fig. 8). The folds cause weak points that can lead to web breaks. The cause of the folds was not known, but the cross-sections suggested that they were created by in-plane compressive failure. The problem was important to many mills and guidance on how they might be avoided was not obvious.

Fig. 8. Crepe wrinkles are the result of compressive failure of paper in the machine-direction, as shown by cross-sections (McDonald 2014)

It was pictured that the wrinkles were created by inter-layer slippage of paper within rolls, which produced in-plane compressive failure of the paper layers. Slippage between layers in the rolls is prevented by friction that depends on the coefficient of the friction of paper (McDonald 1999) and radial compressive stress. This stress is produced by the cumulative tension of overlying layers during winding. Slippage occurs when an applied force exceeds the friction force. The source of this force is discussed below.

In two-drum winders, paper rolls are driven and rest on the winder drums. Force is exerted though the nips between paper and the hard roll. The viscoelastic property of paper causes force in the nip to skew forward, which creates a rolling friction that must be overcome by torque. This torque tends to unwind the roll. If it exceeds layer-to-layer friction, slippage occurs within the paper roll (McDonald and Ménard 1999), which leads to the in-plane compressive failure of paper (McDonald 2014).

This problem was modeled by a force balance to predict the minimum paper-to-paper coefficient friction required to prevent layer-to-layer slippage and buckling of a sheet supported on both sides by an elastic foundation (McDonald 2014). The resulting model gave the minimum required coefficient of friction (μ) as a function of paper roll diameter (R), roll density (ρ), radial compressive stress (σr), the compressive modulus of the roll (Er), and the radius of the winder drums (RD) (Eq. 8):


Equation 8 was simplified by employing previous work that showed that the modulus of a paper roll (ER) exposed to a radial load is proportional to the radial compressive stress (σr) (McDonald et al. 2005)Using the proportionality (Er σr) and for simplicity, an infinite winder drum radius, the minimum coefficient of friction to avoid slippage is given by Eq. 9:


This simple equation provides guidance to control crepe wrinkles and has been found consistent with commercial experience (McDonald 2014).


  1. These examples have illustrated the development and application of several mathematical models in papermaking.
  2. The models addressed a clear problem or need and were based on a conceptual picture of the key variables and how they interact. Formulating the equations brought a deeper understanding, and the results brought new predictive capabilities.
  3. Although there are differing ways of conducting scientific research, this approach has been of great value in the long history of science. Future researchers in pulp and paper science are encouraged to employ it.


Celzard, A., Fierro, V., and Kerekes, R. (2009). “Flocculation of cellulose fibres: New comparison of crowding factor with percolation and effective medium theories,” Cellulose 16, 983-987. DOI: 10.1007/s10570-009-9314-0

Gooding, R. W., and Kerekes, R. J. (1989). “Derivation of performance equations for solid-solid screens,” The Canadian Journal of Chemical Engineering 67(5), 801-805. DOI: 10.1002/cjce.5450670511

Gooding, R. W., and Kerekes, R. J. (1992). “Consistency changes caused by pulp screening,” TAPPI Journal 75, 109-118.

Green, S. I., Zhao, R. H., and Kerekes, R. J. (1998). “Pressure distribution between forming fabrics in blade gap formers: Blades of finite width and fabrics of finite stiffness,” Journal of Pulp and Paper Science 24(2), 60-67.

Hamel, J., Browne, T. C., and McDonald, J. D. (1999). “Calender bubbling: Causes and cures,” Pulp & Paper Canada 100(3), T80.

Kerekes, R. J., Soszynski, R. M., and Tam Doo, P. A. (1985). “The flocculation of pulp fibres,” in: Transactions of the Eighth Fundamental Research Symposium, Oxford, England, pp. 265-310.

Kerekes, R. J., and McDonald, J. D. (1991). “A decreasing permeability model of wet pressing: Theory,” TAPPI Journal 74(12), 150-156.

Kerekes, R. J., and Schell, C. J. (1992). “Characterization of fibre flocculation regimes by a crowding factor,” Journal of Pulp and Paper Science 18(1), 32-38.

Kerekes, R. J., McDonald, E. M., and McDonald, J. D. (2013). “Decreasing permeability model of wet pressing: Extension to equilibrium conditions,” Journal of Science & Technology for Forest Products and Processes 3(2), 46-51.

Kerekes, R. J., and McDonald, J. D. (2020). “Equilibrium moisture content in wet pressing of paper,” TAPPI Journal 19(7), 333-340.

Kubat, J., and Steenberg, B. (1955). “Screening at low particle concentrations,” Svensk Papperstid 58(9), 319-324.

Martinez, D. M., Buckley, K., Jivan, S., Lindstrom, A., Thiruvengadaswamy, R., Olson, J. A., Ruth, T. J., and Kerekes, R. J. (2001). “Characterizing the mobility of papermaking fibres during sedimentation,” in: Transactions of the 12th Fundamental Research Symposium, Oxford, England, pp. 225-254.

Mason, S. G. (1950). “The motion of fibres in flowing liquids,” Pulp and Paper Magazine of Canada 51(10), 93-100.

McDonald, J. D., and Kerekes, R. J. (1991). “A decreasing permeability model of wet pressing: Applications,” TAPPI Journal 74(12), 142-149.

McDonald, J. D., and Kerekes, R. J. (1995). “A decreasing-permeability model of wet pressing with rewetting,” TAPPI Journal 78(11), 107-111.

McDonald, J. D., and Ménard, A. (1999). “Layer-to-layer slippage within paper rolls during winding,” Journal of Pulp and Paper Science 25(4), 148-153.

McDonald, J. D. (1999). “Understanding friction,” in: Proceedings of the Fifth International Conference on Web Handling, Stillwater, OK, USA, pp. 195-215.

McDonald, J. D., Hamel, J., and Kerekes, R. J. (2000). “Design equation for paper machine press sections,” Journal of Pulp and Paper Science 26(11), 401-406.

McDonald, J. D., Hamel, J., and Menard, A. (2005). “Out-of-round paper rolls,” Pulp & Paper Canada 106(5), 35-40.

McDonald, J. D. (2014). “Crepe wrinkle formation during reeling and winding,” Journal of Science & Technology for Forest Products and Processes 4(1), 10-16.

McDonald, J. D., McDonald, E. M., and Kerekes, R. J. (2013). “Impact of press felt design on paper machine press dewatering,” Journal of Science & Technology for Forest Products and Processes 3(2), 52-57.

McDonald, J. D., and Kerekes, R. J. (2017a). “Pragmatic mathematical models of wet pressing in papermaking,” BioResources 12(4), 9520-9537. DOI: 10.15376/biores.12.4.McDonald

McDonald, J. D., and Kerekes, R. J. (2017b). “Estimating limits of wet pressing on paper machines,” TAPPI Journal 16(2), 81-87. DOI: 10.32964/TJ16.2.81

McDonald, J. D., and Kerekes, R. J. (2018). “Rewet in wet pressing of paper,” TAPPI Journal 17(9), 479-487. DOI: 10.32964/TJ17.09.479

Meyer, R., and Wahren, D. (1964). “On the elastic properties of three-dimensional fibre networks,” Svensk Papperstidning 57(10), 432-436.

Nelson, G. L. (1981). “The screening quotient: A better index for screening performance,” TAPPI Journal 64(5), 133-134.

Soszynski, R., and Kerekes, R. J. (1988). “Elastic interlocking of nylon fibers suspended in liquid part 1. Nature of cohesion among fibers,” Nordic Pulp & Paper Research Journal 3(4), 172-179. DOI: 10.3183/npprj-1988-03-04-p172-179

Wigner, E. P. (1960). “The unreasonable effectiveness of mathematics in the natural sciences,” Communications on Pure and Applied Mathematics 13(1), 1-14. DOI: 10.1002/cpa.3160130102

Zhao, R. H., and Kerekes, R. J. (1995). “Pressure distribution between forming fabrics in blade gap formers: Thin blades,” Journal of Pulp and Paper Science 21(3), 97-103.

Article submitted: May 23, 2020; Peer review completed: July 19, 2020; Revised version received and accepted: July 23, 2020; Published: August 5, 2020.

DOI: 10.15376/biores.15.4.7319-7329