Upcycling bleached old corrugated container fibers into creped recycled tissue: Softness–strength trade-offs and resource-substitution potential

Abstract

Limited high-quality bleached recycled fiber (BRF) in North America has increased old corrugated container (OCC) fiber use as a partial substitute. This study examined the effects of replacing BRF with OCC (0 to 100%) on fiber structure, tensile strength, softness, and economics in creped handsheets. As OCC content rose, fiber length (1.084 to 1.710 mm) and coarseness (107 to 186 mg·km⁻¹) increased, while kink index decreased (1.961 to 1.501). These changes improved bulk and reduced softness parameters TS7 and TS750. Tukey HSD identified four statistical groups for TS7 and three for TS750 (α = 0.05). The tensile index remained stable up to 75% OCC but dropped 31.9% at 100%, with a 35.0% increase in bulk. Sensory tests ranked 100% OCC as softest, aligning with instrumental data. A Monte Carlo simulation showed 77.7% of cases had lower costs with treated OCC, median savings of USD 80/odt, and a break-even price of USD 252/odt. Factors affecting competitiveness were yield and delivered price. Substituting 25 to 75% OCC reduces costs with minimal strength loss; total replacement maximizes savings and bulk but reduces tensile strength.

Full Article

Upcycling Bleached Old Corrugated Container Fibers into Creped Recycled Tissue: Softness–Strength Trade-Offs and Resource-Substitution Potential

Fernando Urdaneta ,^a Ronald Marquez ,^b Anibal Moran ,^a Isabel Enriquez-Medina ,^a Dipta Mozumder ,^a Isabel Urdaneta ,^a Joel J. Pawlak ,^a Hasan Jameel ,^a Daniel Saloni ,^a and Ronalds W. Gonzalez ^a,*

Limited high-quality bleached recycled fiber (BRF) in North America has increased old corrugated container (OCC) fiber use as a partial substitute. This study examined the effects of replacing BRF with OCC (0 to 100%) on fiber structure, tensile strength, softness, and economics in creped handsheets. As OCC content rose, fiber length (1.084 to 1.710 mm) and coarseness (107 to 186 mg·km⁻¹) increased, while kink index decreased (1.961 to 1.501). These changes improved bulk and reduced softness parameters TS7 and TS750. Tukey HSD identified four statistical groups for TS7 and three for TS750 (α = 0.05). The tensile index remained stable up to 75% OCC but dropped 31.9% at 100%, with a 35.0% increase in bulk. Sensory tests ranked 100% OCC as softest, aligning with instrumental data. A Monte Carlo simulation showed 77.7% of cases had lower costs with treated OCC, median savings of USD 80/odt, and a break-even price of USD 252/odt. Factors affecting competitiveness were yield and delivered price. Substituting 25 to 75% OCC reduces costs with minimal strength loss; total replacement maximizes savings and bulk but reduces tensile strength.

DOI: 10.15376/biores.21.4.11129-11166

Keywords: Upcycling; Softness; Sustainability; Fibers; Tissue; Recycling; Papermaking

Contact information: a: Department of Forest Biomaterials, College of Natural Resources, North Carolina State University, Campus Box 8005, Raleigh, NC 27695‑8005 USA; b: LEPAMAP-PRODIS Research Group, University of Girona. C/ Maria Aurèlia Capmany, 61 – 17003 Girona (Spain);

* Corresponding author: [email protected]

INTRODUCTION

The global tissue and towel industry relies heavily on recycled fibers to meet growing consumer demand while reducing dependence on virgin wood resources. Recycled fiber currently accounts for a substantial share of total fiber input in tissue manufacturing worldwide, especially for away-from-home products, and North America is among the most recycled-fiber-dependent tissue markets globally, with most manufacturers incorporating recovered paper into their furnish (Moore and Cornish 2020). For tissue grades requiring acceptable brightness, softness, and cleanliness, the preferred recycled furnish has traditionally been sorted office paper (SOP) and other deinked high grades derived from office and copy papers, particularly uncoated freesheet (UFS) (Christopher 2020; Moore and Cornish 2020). After deinking, these bleached recycled fibers (BRF) provide the combination of fine fiber dimensions, low lignin content, and adequate swelling capacity needed to achieve acceptable tactile and mechanical properties in a recycled tissue furnish (Moore and Cornish 2020; U.S. Environmental Protection Agency 2000). Tissue softness is a multidimensional property encompassing surface smoothness and bulk compliance, both of which are closely tied to fiber morphology and network structure (Wang et al. 2018; de Assis et al. 2019; Pawlak et al. 2022; Urdaneta, I., et al. 2025). High-quality deinked bleached fibers have historically supplied the fine, flexible fiber fraction that underpins softness performance in recycled tissue grades (Hollmark and Ampulski 2004; D’Olivo 2022; Pawlak et al. 2022).

However, the sustained availability of these high-quality deinked fiber streams faces growing structural pressure from the global shift toward digital communication. The widespread adoption of remote work, e-billing, electronic publishing, and digital office environments has led to persistent, accelerating declines in the production and consumption of printing and writing papers, particularly UFS (Berg and Lingqvist 2019; Moore and Cornish 2020). Globally, printing and writing paper consumption has declined while tissue consumption has grown over the past two decades, widening the gap between recoverable high-grade fiber supply and tissue industry demand (Berg and Lingqvist 2019; Moore and Cornish 2020). In North America, this structural contraction has continued, and trade-press and association reports highlight ongoing year-over-year decreases in printing and writing paper shipments (Moore and Cornish 2020; TissueOnline 2023; World Paper Mill 2026). The COVID-19 pandemic accelerated this trend, as extended remote work has mandated sharply reduced office paper generation and SOP collection volumes, forcing tissue mills to either accept lower-quality fiber alternatives or pay elevated premiums for the diminishing supply of preferred deinked grades (Christopher 2020; Moore and Cornish 2020; TissueOnline 2023).

Given these supply constraints, Old Corrugated Containers (OCC) are the most abundant and logistically accessible recovered-fiber category in North America. OCC is the most widely recycled paper grade globally, with repulping yields on the order of 80%. Domestic mill demand for OCC is projected to continue growing alongside expanding recycled containerboard and packaging capacity (Zambrano et al. 2021; Quinn 2024). The composition of containerboard products, however, varies across regions. North America is unique among major producing regions in maintaining a comparatively larger share of virgin and mixed-furnish products. While products classified as virgin often contain a proportion of recycled fiber for economic and performance benefits, they are generally categorized as mixed-furnish grades. In contrast, containerboard production in most other regions is predominantly based on recycled fiber (Veijalainen 2022).

These regional differences in furnish composition can influence the characteristics of the OCC stream and, consequently, its suitability for tissue manufacturing. Variations in fiber origin and recycling history affect fiber morphology and chemical composition, all of which play a critical role in determining the performance of recovered fibers in downstream applications. Unlike high-grade deinked pulps, OCC fibers have a fundamentally different structural profile that poses well-documented challenges for tissue applications. The relatively high residual lignin content and extensive recycling history of OCC make its fibers inherently stiff and prone to hornification, which reduces fiber swelling capacity, bonding potential, and flexibility with each successive recycling cycle, characteristics broadly regarded as detrimental to tissue softness and tensile strength (Zambrano et al. 2021, 2022). Recent investigations of upcycled OCC-based brown tissue demonstrate that acceptable absorbency and softness can be engineered by modifying furnish design and processing, but these products still differ markedly from conventional bleached recycled tissue in fiber morphology and network structure (Zambrano et al. 2021, 2022). Prior work within the Sustainable and Alternative Fiber Initiative (SAFI) at North Carolina State University has shown, through comparative evaluations of wood and non-wood market pulps, that furnish composition strongly governs bulk, softness, and tensile performance in creped tissue, with fiber morphology and bonding capacity as primary mechanistic levers (de Assis et al. 2018, 2019; Urdaneta, F., et al. 2024, 2025; Urdaneta, I., et al. 2025).

Despite OCC’s abundance and logistical accessibility, quantitative data on how progressively substituting BRF with OCC affects instrumental and sensory softness under creped conditions remain limited, particularly when basis weight is controlled, and bulk, TS7, TS750, and tensile properties are measured on the same handsheets. The Tissue Softness Analyzer (TSA) provides two widely used softness indices: TS7, which reflects surface softness and correlates inversely with surface roughness, and TS750, which quantifies bulk stiffness of the sheet structure (EMTEC Electronic GmbH 2025; Urdaneta, I., et al. 2025). A validated study confirmed that TS7 closely correlates with human panel scores for bath tissue, establishing it as a reliable proxy for perceived softness (Wang et al. 2018). More recently, a statistical model combining TSA parameters with fiber morphological data achieved high correlation with panel handfeel in tissue containing alternative fibers, highlighting the value of integrated instrumental–sensory approaches to furnish design (Urdaneta, F., et al. 2024, 2025; Urdaneta, I., et al. 2025). Despite the demonstrated utility of this framework, it has not been applied to creped laboratory handsheets made from progressive OCC and bleached recycled fiber blends, the furnish design space most relevant to mills navigating declining SOP availability.

This study had three objectives: first, it quantified how replacing BRF with chelation–peroxide-treated OCC at ratios from 0 to 100% affected fiber morphology, volume, tensile properties, TSA response, and softness of handsheets. Second, it assessed the consistency between instrumental softness measures and sensory rankings, with interpretation constrained by the use of five furnish-level means. Third, it developed a furnish-level techno-economic assessment, including uncertainty in recovered-paper prices, yields, chemical needs, utilities, allowances, and disposal costs via Monte Carlo simulation. Based on the greater stiffness, coarseness, and hornification of OCC fibers, increasing OCC content was initially expected to reduce tissue softness. However, preliminary TSA measurements showed an unexpected decrease in TS7 with increasing OCC content. This counterintuitive response motivated the sensory panel evaluation to determine whether the instrumental trend was also reflected in perceived softness rather than being limited to the TSA measurement. The sensory and instrumental results subsequently showed directional agreement in the softness response. The economic analysis aimed to find cost and break-even conditions, not profitability. Together, the experiments and economic model provide a basis for identifying promising OCC substitution ranges for further validation.

EXPERIMENTAL

Materials

Bleached recycled fiber (BRF) was obtained as wet pulp from a North Carolina tissue mill that primarily processes deinked recycled furnish from sorted office paper (SOP) and uncoated freesheet (UFS). The mill’s deinking process includes centrifugal cleaners, pressure screens, and a high-consistency disk disperser, which applies mechanical shear to reduce residual ink, stickies, and wax particles to sub-micron sizes. The pulp then undergoes oxidative bleaching to enhance brightness and fiber surface hydrophilicity. The BRF sample, collected at 6% solids after the disperser and bleaching tower but before stock refiners, reflects a fully cleaned and bleached but mechanically unrefined state, preserving fiber morphology and surface properties. The pulp was centrifuged to approximately 30% consistency, stored at 2 °C to inhibit microbial activity, then warmed and gently disintegrated in deionized water for stock preparation.

The OCC material was sourced from corrugated boxes in the Northwestern US, manually cleaned to remove tapes, labels, and debris, then pulped in deionized water, screened through a 0.15 mm slotted screen, and bleached via a two-stage chelation–peroxide process. The chelation stage used EDTA (1%) at 50 °C for 90 min to remove transition-metal ions, followed by a peroxide stage with 4% H₂O₂, 1.5% NaOH, and 1.6% sodium silicate at 90 °C for 120 min. Chemical products used in this work were EDTA (ethylenediaminetetraacetic acid, 99.5% ACS reagent; Acros Organics), H₂O₂ (34%, FCC; Lab Alley), NaOH (50% w/w certified solution; Fisher Chemical), and sodium silicate supplied as anhydrous sodium metasilicate (technical grade; Thermo Scientific). No refining or additional mechanical treatment was applied. The OCC had a freeness of 520 mL CSF, while the BRF had a freeness of 435 mL CSF. A complete oven-dry mass balance across screening, bleaching, washing, and reject streams was not determined; therefore, the bleaching treatment is reported as a processing condition rather than as a measured yield analysis.

Stock Preparation

For stock preparation, the 30% consistency BRF was diluted with deionized water and disintegrated to a working stock consistency of approximately 3% to ensure full fiber wetting and dispersion. The screened OCC pulp was similarly diluted to 3% consistency. BRF and OCC stocks were then combined gravimetrically, on an oven-dry basis, to achieve the target blend ratios for Samples A through E. All furnishes were adjusted to a forming consistency of 0.3% with deionized water and gently stirred to ensure homogeneous fiber distribution immediately prior to sheet formation. No refining, beating, or chemical additives were used in stock preparation.

Fiber Morphology Analysis

Fiber morphological characterization was carried out using a Fiber Quality Analyzer (FQA; OpTest Equipment Inc., Hawkesbury, ON, Canada) in flow-through mode. Before measurement, each pulp sample was diluted to approximately 0.02% (w/v) in deionized water to ensure single-fiber flow and reduce fiber overlap during optical detection. Three replicate measurements were taken for each furnish sample. The following morphological parameters were recorded for each replicate: length-weighted mean fiber length (L_w, mm), arithmetic mean fiber width (µm), ISO coarseness (mg/km), fiber wall thickness (FWT, µm), arithmetic mean curl index, arithmetic mean kink index, and length-weighted percent fines (%). Length and fines were measured within the 0.1 to 5.0 mm range, with fines defined as fibers shorter than 0.2 mm (ISO cutoff). Curl and kink indices were calculated for fibers between 0.5 and 5.0 mm. Coarseness was measured following ISO/FDIS 23713:2005 (E).

Handsheet Formation and Conditioning

Laboratory handsheets were formed in a British Sheet Machine manufactured by Robert Mitchell Inc. (Canada) and distributed by Testing Machines Inc. (Amityville, NY, USA). The target basis weight of the uncreped handsheets was 30 g/m² for all furnishes. After formation, each wet handsheet was couched under two blotter papers, sandwiched between them, sealed in polyethylene bags, and stored for 24 h to equilibrate moisture content across the batch prior to creping. For each furnish composition, a single stock batch was prepared and used for handsheet production.

Laboratory-Scale Creping and Post-Creping Stretching

Creping was performed using a custom-built high-speed laboratory creping simulator developed for tissue research and described in detail by Hämäläinen et al. (2016). The setup included a horizontally oriented, electrically heated sled representing the Yankee dryer surface, a pressurized adhesive spray system, and a stationary creping blade positioned downstream along the sled’s path. The sled was heated to 105 °C, and the creping speed was maintained at 600 m/min, consistent with conditions from de Assis et al. (2019). A PVA-based adhesive solution was uniformly sprayed onto the sled at a controlled rate to facilitate bonding with the wet web. An uncreped handsheet was transferred onto the moving sled using a press roll, ensuring close contact with the heated, adhesive-coated surface. As the sheet moved into the creping zone, it was impacted by a stationary blade set at an 80° angle relative to the sled’s direction, causing interfacial delamination, buckling, and microfold formation (Hämäläinen et al. 2016). The creped sheet was immediately collected and dried flat to prevent geometric distortion. Between tests, the sled and blade areas were cleaned to maintain consistent creping conditions. After creping, each handsheet underwent a controlled stretching process to relieve residual stresses and promote uniform micro-fold development (Vieira et al. 2022). This procedure was performed on the same testing machine used for tensile tests, following TAPPI T 581 om-22 (2022). The sheets were cut to fit the machine’s grip width, averaging 7.4 cm across all samples. Prior to stretching, the gauge length at failure was determined from preliminary tests; the stretching gap was set to 90% of this length to induce sub-failure deformation. Each specimen was held at this extension for 20 min before release. Finally, all samples were conditioned at 23 °C and 50% RH as per TAPPI T 402 sp-21 (2021) before measuring their properties. These operating conditions were kept constant for all furnishes to minimize sample-to-sample differences in creping. Nevertheless, small variations in sheet adhesion, transfer, and microfold formation may contribute to the observed within-furnish variability.

Physical and Mechanical Testing

All property measurements were conducted on conditioned (at 23 °C and 50% relative humidity) handsheets (n = 5 replicates per furnish). Basis weight of the creped sheets was determined following TAPPI T 410 om-23 (2023) (except for the number of samples, which in this study was 5) using accurately cut test pieces of known area. Caliper (thickness) was measured according to ISO 12625-3 (2014), and bulk was calculated as the ratio of caliper to basis weight. Because the handsheets were formed isotropically in the British-style sheet former, they possess no defined machine or cross-machine direction; tensile properties were therefore measured in the creping direction on a universal testing machine following TAPPI T 581 om-22 (2022). Dry tensile properties were measured in the creping direction using a horizontal tensile tester (Testing Machines Inc., New Castle, DE, USA; model 84 56 00 0003, software V1.0.1, 1000 N load range). Rectangular strips (50 mm width, 100 mm gauge length) were tested at 25 mm/min crosshead speed with a 5 N force threshold, 0 N pre load, and a 95% break sensitivity. For each furnish, the instrument calculations of tensile strength, tensile index, stretch, tensile energy absorption (TEA), stiffness, elastic modulus, and breaking length used the measured grammage and thickness of that sheet, as reported in the physical properties section.

The effective crepe ratio (CR) for each furnish was calculated from the change in basis weight before and after creping, as follows,

(1)

where BW_uncreped is the basis weight of the uncreped handsheet (target 30 g/m²) and BW_creped is the basis weight of the creped sheet. This approach is consistent with the definition of percent crepe used in tissue manufacturing, where creping shortens the sheet and concentrates fiber mass per unit area, causing the creped basis weight to exceed the uncreped target (de Assis et al. 2019; Padley 2022).

Brightness and color were measured using a Technidyne ColorTouch X reflectance spectrophotometer (Technidyne Corporation, New Albany, IN, USA) with a D65 illuminant and 10° standard observer setup (D65/10°). The measurements included ISO brightness, CIE L*, a*, and b* color coordinates, the CIE whiteness index (WI), tint, and yellowness index (YI; ASTM E313).

Softness Measurements

Instrumental softness was evaluated using a Tissue Softness Analyzer (TSA; emtec Electronic GmbH, Leipzig, Germany) following the approach described by Prinz et al. (2021), Wang et al. (2018), and Prinz et al. (2021), and validated against human panel scores by Wang et al. (2019). Five replicates per furnish were tested with the sled-contact surface facing the TSA blades. For each measurement, the TSA recorded the full acoustic response spectrum over a frequency range of 0 to 23,950 Hz, expressed in decibels (dB). TS7, which reflects high-frequency surface vibration associated with fiber-level roughness and has been shown to correlate linearly with perceived softness in bath tissue (Wang et al. 2018), was used as the primary softness indicator. TS750, which quantifies the low-frequency structural response related to bulk stiffness of the sheet, was recorded as a secondary parameter (Prinz et al. 2021). Sample preparation, clamping, and testing conditions followed the instrument manufacturer’s recommendations (EMTEC Electronic GmbH 2025).

Sensory Panel Testing

A forced-choice sensory panel was conducted to assess perceived softness differences among the five furnishes (Samples A through E) using a round-robin pairwise comparison design (Urdaneta, I., et al. 2025). Twenty-six panelists participated in the evaluation. Each panelist assessed all 10 pairwise combinations of the five samples (A vs. B, A vs. C, A vs. D, A vs. E, B vs. C, B vs. D, B vs. E, C vs. E, D vs. E, and C vs. D), resulting in 260 total comparisons. Panelists were presented with blinded, randomized pairs of creped tissue samples and asked to identify the softer sample in each pair by handfeel. All samples were conditioned at 23 °C and 50% relative humidity before evaluation, and sessions were conducted in a controlled laboratory environment.

Paper Properties Data Analysis

Descriptive statistics (mean and standard deviation) were computed for all physical, mechanical, and softness properties. To provide an integrated index of the overall TSA spectral response of each furnish, the area under the full TSA spectrum (AUC) was computed by trapezoidal numerical integration of the mean TSA response (dB) over the complete measured frequency range from 0 to 23,950 Hz, according to Eq. 2 (a fully worked numerical example is provided in Supplementary Information, Section S1):

(2)

A(f) is the mean TSA response in dB at frequency f (Hz), evaluated over the recorded frequency range. For each furnish, the replicate TSA spectra were first averaged point-by-point, and the area under the resulting mean spectrum (AUC), in dB·Hz, was then calculated by trapezoidal numerical integration. The AUC was used as a descriptive furnish-level index of the integrated TSA spectral response and was not treated as an independent replicate-level response for inferential statistical testing. Because the TSA signal is reported in decibels, the AUC should not be interpreted as physical acoustic energy. Numerical integrations were performed in Python 3.10 using numpy.trapezoid on the mean dB array. For sensory evaluation, the full 26-panelist dataset was treated as the primary analysis. A secondary sensitivity analysis excluded panelists showing systematic disagreement with the aggregate ordering using Kendall’s τ ≤ −0.2 as a pragmatic screening criterion; the rationale and calculation are provided in Supplementary Information, Section S7.2. All data processing workflows were executed utilizing the pandas library.

One-way ANOVA was used to test differences among five furnishes for physical, mechanical, and softness properties. Before each ANOVA, Levene’s test assessed homogeneity of variance. Levene’s test was non-significant (p ≥ 0.05), and pairwise comparisons used Tukey’s HSD. CLD (compact letter display) diagrams from the sweep algorithm at α=0.05 indicate that means sharing a letter are not different. All statistical modeling was performed in Python 3.10 using scipy.stats and numpy libraries.

All Pearson correlations were calculated across the five furnish-level means representing a gradient. Since these points are sequential substitutions, the correlation coefficients describe gradient co-variation and should be seen as indicative, not confirmatory. Properties such as TS7, tensile index, and perceived softness, which respond non-monotonically to OCC substitution, are especially prone to inflated r values, so their associations should be interpreted cautiously.

Three complementary statistical methods were applied to the sensory panel data. A Bradley–Terry paired-comparison model was fitted iteratively to the full win–loss matrix to estimate a latent softness strength parameter for each furnish (Bradley and Terry 1952). This probabilistic model assumes that the probability of furnish i being preferred over furnish j is a function of their respective latent strength scores and converges to maximum likelihood estimates through iterative scaling (Bradley and Terry 1952). The resulting latent scores provide a coherent probabilistic ranking of the five furnishes on a continuous scale.

Thurstone Case V scaling was used to convert the observed pairwise win proportions into a continuous psychological interval scale (Thurstone 1927). Pairwise proportions were transformed using the inverse normal (z) function under Case V assumptions of equal and independent discriminable dispersions, and scale values were obtained as column means of the resulting z-score matrix (Thurstone 1927). The scale values quantify the magnitude of perceived softness distances between furnishes and are directly interpretable as interval-level measurements.

Two-tailed binomial sign tests (H₀: p = 0.50, α = 0.05) were applied to each of the 10 pairwise comparisons to assess whether the observed preference for one furnish over another could be attributed to chance (Siegel and Castellan 1988). For each pair, the number of votes for the majority “softer” sample was compared to the binomial distribution with n equal to the number of valid panelist responses.

Additional supplementary analyses were performed to assess panel behavior and robustness and are reported in full in the Supplementary Information. These include Kendall’s coefficient of concordance (W) to quantify overall agreement among panelists (Kendall, 1975) (Supplementary Information, Section S7.3), intransitive triad analysis to evaluate individual-level consistency in paired comparisons (Section S7.1), and post-hoc power calculations to evaluate the adequacy of the panel size at the observed effect sizes (Section S6). Full formulas, the complete win matrix, intermediate z-score matrices, convergence traces, and numerical outputs for all analyses are provided in the Supplementary Information (Sections S1–S7).

Screening-Level Furnish-Cost Model: Market and Process Assumptions

The recovered-paper process assumed an input cost of USD 250 per odt with a 70% yield, modeled as a scenario rather than a typical case (Zambrano et al. 2021). The delivered paper price followed a triangular distribution from USD 200 to 325 per odt, with a most probable value of USD 250 (Morris 2018; Zambrano et al. 2021). Usable yield at the headbox was represented by a distribution from 60% to 85%, most likely 70%. Due to limited data, costs for repulping, cleaning, deinking, and wastewater were estimated with a triangular distribution ranging from USD 40 to 140 per odt, with a mode at USD 80, while fixed costs for labor, maintenance, and capital ranged from USD 20 to 100 per odt, most likely USD 50. For the treated-OCC route, the delivered price ranged from USD 100 to 300 per odt, with a mode at USD 175 (Workman 2024; Zambrano et al. 2021), and usable yield from 65% to 90%, most likely 75%, accounting for processing losses. All values are nominal USD, and fiber measures are based on oven-dry basis. Dry-reject disposal costs were assigned route-specific triangular distributions: USD 30 to 80 per odt reject (mode USD 55) for OCC screening rejects, reflecting standard landfill tipping fees for dry, low-ash residuals (Tanoh et al. 2022); and USD 45 to 150 per odt reject (mode USD 85) for deinking sludge, reflecting the increased dewatering and handling costs documented for sludge disposal (Abdullah et al. 2015).

Goal, Functional Unit, and Comparison Basis

Two recovered-fiber routes ending at the same functional unit were assumed: one using OCC, and the other involving the purchase and processing of higher-grade recovered paper. The system boundary included recovery, yield losses, chemicals, energy, water, effluent handling, reject disposal, and labor costs, but excluded drying, converting, product performance, and transportation costs. The goal was to determine furnish-level cost ranges and uncertainties, not to provide a complete mill-cost or profit estimate.

Economic Model

For each route j, the cost of one oven-dry tonne of usable headbox-ready fiber was calculated as,

(3)

where C_j is the cost of usable headbox-ready fiber for route j in USD/odt usable fiber; P_j is the delivered recovered-paper purchase price in USD/odt input; C_conversion,j is the route-specific chemical, energy, water, and effluent cost in USD/odt input; C_fixed,j is a screening allowance for labor, maintenance, and annualized capital in USD/odt input; Y_j is the overall usable dry-mass yield expressed as a fraction; and CD is the dry-reject disposal cost in USD/odt reject. Division by Y_j expresses input-basis costs per tonne of usable output, while the final term accounts for the disposal of dry material rejected during processing. Reject disposal cost, CD, was modeled separately for the two routes (CD,OCC and CD,RCP) because the reject streams differ physically: OCC-route rejects are coarse, low-moisture screening residuals (staples, tape, film, shives), whereas conventional-route rejects are deinking sludge that must be dewatered before disposal.

Costs for intermediate furnish blends (25, 50, and 75% treated OCC) were determined by mass-weighting the underlying route-level yield, delivered price, and conversion costs according to the OCC mass fraction x, then solving Eq. 3 individually for each blended route. A more straightforward approach involves linear interpolation between the two endpoint route costs as shown in Eq. 4.

(4)

The economic advantage of the treated-OCC route was calculated using Eq. 5:

(5)

Positive values show that the treated-OCC route had a lower modeled furnish cost, while negative values indicate the conventional higher-grade recovered-paper route was cheaper. This difference represents the modeled furnish-level cost advantage available to absorb additional costs not considered in the present screening-level model. It does not represent operating profit.

Treatment-Chemical and Utility Basis

In the Q–P-treated OCC process, the active chemical cost was calculated by multiplying each chemical’s load by its unit price separately, rather than using a combined factor. The typical active loads per odt of OCC were 10 kg EDTA, 40 kg H2O2, 15 kg NaOH, and 16 kg sodium silicate. In the base Monte Carlo model, these loads and prices were varied independently. Theoretical purchased heat was estimated based on slurry consistency, a total process-temperature rise of 90 K, heating efficiency, and heat recovery. The model did not fully include commercial reagent concentrations, solution water, storage, freight, chemical losses, pumping, agitation, screening electricity, dewatering electricity, or supplier-specific delivery conditions. These should be incorporated when specific mill or supplier data are available.

Monte Carlo Uncertainty Propagation and Global Sensitivity Analysis

Uncertainty was propagated using 100,000 Monte Carlo trials with a fixed seed (20260730). Triangular distributions were used because inputs were defined by plausible minimum, mode, and maximum values rather than large empirical datasets suitable for parametric fitting. Each trial sampled one value per uncertain price, yield, chemical load/price, utility, fixed-cost, and reject-disposal input, then computed route-specific usable-fiber costs, the economic advantage of treated OCC, and the break-even OCC purchase price.

Inputs were sampled independently because the screening dataset lacked paired observations to justify a dependence structure. As a result, the probabilities are conditional on the specified marginals and the assumption of independence, not actual market frequencies. In real-world applications, the recovered-paper price might correlate with quality and yield, and treatment intensity could be linked to chemical use, energy demand, yield, brightness, and effluent costs. Once paired mill data are available, a scenario with correlated inputs should be considered. Reject disposal costs were also sampled independently for each route due to their unique waste characteristics and disposal methods.

Reported outputs include the mean, median, standard deviation, P10, P90, and probability that treated OCC costs less than the conventional route. The break-even OCC price was found in each trial by solving Eq. (3) for the OCC input price equating both route costs, holding all other sampled parameters fixed; a negative value indicates treated OCC would be non-competitive even at zero purchase cost under that trial’s assumptions.

Global sensitivity used Spearman rank correlations between each input and the simulated economic advantage, chosen for capturing monotonic (not necessarily linear or normal) relationships. Coefficients indicate influential assumptions within the specified ranges, not causal effects, and rankings are conditional on those ranges and the independence assumption. Convergence should be checked by running estimates of the median advantage, P10/P90, win probability, and median break-even price.

RESULTS AND DISCUSSION

Fiber Morphology

The fiber quality attributes of the five BRF/OCC furnish types are shown in Table 1 and in Fig. 1 (as normalized values to the 100% BRF benchmark).

Table 1. Fiber Quality Analyzer (FQA) Morphological Properties of the Five BRF/OCC Furnishes

Consistent with the findings of Reynolds et al. (2025), the morphological analysis indicated that the OCC pulp contained a substantial proportion of softwood fibers, likely originating from the linerboard component of containerboard manufactured in the United States (Reynolds et al. 2025; Veijalainen 2022). As OCC substitution increased, all measured morphological parameters changed systematically, reflecting the fundamental structural differences between the hardwood-dominated BRF furnish and the softwood-enriched OCC pulp.

The length-weighted mean fiber length (Lw) steadily increased from 1.08 mm in Sample A (100% BRF) to 1.18, 1.33, 1.49, and 1.71 mm in Samples B through E. This trend reflects the gradual addition of longer-fibered OCC fractions into the furnish. Notably, Sample E (100% OCC) had an L_w nearly 50% higher than Sample A, aligning with the dominance of softwood kraft fibers in U.S. corrugated containerboard.

Coarseness increased from 107 mg/km in Sample A to 186 mg/km in Sample E, representing a 1.7-fold rise across the blend gradient. This shift reflects the substitution of finer, lower-coarseness deinked hardwood fibers—typical of SOP/UFS-derived BRF—with coarser, higher-cell-wall-mass softwood fibers from OCC. Fiber wall thickness (FWT), estimated using the FQA optical model, showed a similar trend, increasing from 1.389 µm (Sample A) to 2.120 µm (Sample E). These findings support the coarseness results and suggest that the increased cell wall mass per unit length in OCC-rich samples genuinely reflects differences in cell wall structure rather than measurement error.

The length-weighted fines content dropped from 5.39% in Sample A to as low as 3.77% in Sample E, marking a 30.1% overall decrease. The reduction was not strictly consistent, as Sample D (4.59%) showed a slight increase compared to Sample C (4.25%). This minor non-monotonic variation at 75% OCC likely results from natural variability in the fines population at intermediate blend ratios or small differences in the mechanical dispersion history of individual stock batches. However, this does not affect the overall trend of decreasing fines content with increasing OCC levels.

Fig. 1. Fiber-morphology means are normalized to the 100% BRF control. Lw, length-weighted fiber length; FWT, apparent FQA fiber-wall thickness

Both the curl index and kink index decreased steadily as OCC content increased. The curl index dropped from 0.131 in Sample A to 0.101 in Sample E, representing a 22.9% decrease. The kink index fell from 1.961 to 1.501, a 23.5% reduction. This suggests that OCC fibers are, on average, straighter and less kinked than BRF fibers. This aligns with the understanding that hornified, high-coarseness softwood fibers are less flexible and stiffer, making them more resistant to shear deformation during stock preparation.

Tensile Properties of Creped Handsheets

Because fiber length increased and kink index decreased with OCC substitution (Table 1; Fig. 1), tensile properties were examined next to determine whether these morphological shifts affected sheet bonding and strength.

The tensile properties of the five creped furnishes are presented in Table 2. Tensile strength ranged from 0.209 kN/m for Sample E to 0.321 kN/m for Sample B, with Samples A, C, and D occupying an intermediate range of 0.296 to 0.314 kN/m. Tensile index, which normalizes tensile strength by grammage to allow direct comparison across furnishes, followed a broadly similar pattern: Samples A through D exhibited closely grouped values between 7.211 and 7.993 kNm/kg, while Sample E (100% OCC) was markedly lower at 4.912 kNm/kg, representing a reduction of approximately 32% relative to the 100% BRF benchmark (Sample A). The pronounced decrease in tensile index at 100% OCC is consistent with reduced inter-fiber bonding in the OCC-rich network. The OCC furnish contained substantially coarser fibers with greater apparent fiber-wall thickness and fewer fines, which likely reduced fiber conformability and bonded area per unit mass. This is consistent with its higher freeness (520 mL CSF versus 435 mL CSF for BRF). Repeated recycling and hornification may decrease fiber swelling and flexibility, so at full OCC substitution the loss in bonding capacity appears to outweigh the benefit of the longer fibers.

Table 2. Tensile Properties of Creped Handsheets Measured in the Creping Direction (Mean ± SD)

Properties trade-offs are summarized in Table 3 and Fig. 2. Samples A through D did not exhibit a systematic decrease in tensile index with increasing OCC content. Sample C (50% BRF / 50% OCC) recorded the highest mean tensile index in the series at 7.99 kN·m/kg, followed by Samples B (7.77 kN·m/kg), D (7.58 kN·m/kg), and A (7.21 kN·m/kg). Only Sample E (100% OCC) fell substantially below this range, at 4.91 kN·m/kg. Tukey HSD confirmed that Samples A through D were statistically indistinguishable (shared letter a, p = 0.006 overall), while Sample E was assigned to a distinct lower group (letter b). This supports the conclusion that up to 75% OCC substitution did not significantly compromise tensile index under the creping conditions applied. This pattern among Samples A through D does not follow directly from the fiber coarseness and stiffness properties of OCC alone and may instead reflect the combined influence of creping-induced microfold redistribution and residual basis weight variation at intermediate blend ratios.

Table 3. Trade-offs of Replacing BRF with OCC

Fig. 2. Furnish-dependent TS7, TS750, specific volume, and tensile index. Points are means and error bars are ±1 SD (n=5 as reported). Lower TS7 and TS750 correspond to lower instrument-specific surface and structural responses; these parameters are not direct modulus or energy measurements.

Stretch at break was the most variable tensile property across the furnish gradient, ranging from 8.9% for Sample C to 18.4% for Sample B, with no consistent trend across the blend series. Sample B showed significantly higher stretch than all other samples, a pattern confirmed across replicates (coefficient of variation 13.31%), making measurement error unlikely. The higher stretch in Sample B may be due to the specific microfold geometry and residual stress from creping, since stretch in creped tissue is mainly influenced by microfold development rather than fiber extensibility alone. Tensile energy absorption (TEA) was also highest for Sample B (31.64 ± 5.483 J/m²), reflecting the combined effects of adequate tensile strength and high stretch. Samples A, C, D, and E had TEA values in a narrower range (14.4 to 19.7 J/m²), with Sample E the lowest at 15.26 ± 2.657 J/m².

Physical Properties of Creped Handsheets

Given that tensile index was preserved up to 75% OCC despite the morphological changes described above (Table 2), bulk and caliper were evaluated to establish whether sheet architecture, rather than fiber bonding, accounted for the structural transition observed at 100% OCC.

Table 4 displays the physical properties of the five creped furnishes. The basis weight exhibited a significant overall effect (one-way ANOVA, p = 0.019), but the values were confined within a narrow range (39.0 to 42.6 g/m²). Tukey HSD post-hoc tests found no significant pairwise differences between any two furnishes, indicating that variations in grammage were not the main factor influencing softness. The caliper was similar among Samples A to D (223–237 µm; all labeled as “b”) but increased noticeably for Sample E (315.2 ± 10.18 µm; labeled as “a”), revealing a structural change specific to 100% OCC content.

The specific volume followed the same two-tier pattern as caliper. Samples A through D formed a single statistically indistinguishable group (5.47 to 6.08 cm³/g; all “b”), while Sample E (100% OCC) was significantly higher at 7.41 ± 0.25 cm³/g (“a”; one-way ANOVA with Tukey HSD, p < 0.001). Overall, these findings suggest that both bulk and caliper stayed structurally stable across all BRF:OCC blend ratios up to 75% OCC substitution. A notable structural change only occurred at 100% OCC, aligning with the reduced bonding and increased coarseness typical of heavily recycled OCC fibers at full substitution.

The effective crepe ratio (CR) provided additional insight into the structural mechanics of the sheet. An omnibus ANOVA (p = 0.017) revealed a significant overall effect of furnish on CR. However, Tukey HSD tests showed no significant pairwise differences, with all samples grouped under the same letter “a”, despite their numerical range from 23.1% (Sample D) to 29.4% (Sample E). Unlike caliper and bulk, which shifted abruptly only at the 100% OCC level, CR displayed a non-monotonic, U-shaped pattern. Both extremes (100% BRF and 100% OCC) had the highest crepe ratios, while the intermediate blends at 50:50 and 25:75 BRF:OCC had the lowest values.

The stability of bulk and caliper across Samples A through D suggests that softness differences among these four furnishes are unlikely to be explained by gross sheet architecture alone, pointing instead toward fiber-level mechanisms, such as fiber length, coarseness, and fines content, as the more probable drivers of the TS7 and sensory panel trends discussed subsequently.

Table 4. Physical Properties of Creped Handsheets (Mean ± SD, n = 5 per Furnish)

The optical data in Table 5 show a systematic decline in brightness and CIE from sample A to sample E, accompanied by continuous increases in a*, b*, and yellowness index, while the CIE whiteness index transitions from strongly positive to markedly negative values. This trend is consistent with a progressive increase in OCC content, where higher lignin and chromophore contents increase light absorption and shift the shade toward red-yellow, thereby lowering perceived whiteness under D65/10° illumination.

Table 5. Instrumental Brightness and Color Parameters of Samples A–E, Measured with a D65/10° Spectrophotometer

Instrumental Softness (TSA)

Because bulk and caliper remained statistically similar across Samples A to D, the corresponding TSA changes cannot be explained by bulk and caliper alone, suggesting that fiber- and network-level characteristics also contributed to the observed acoustic response. With fiber length, coarseness, and kink index changing systematically with OCC content (Table 1), while tensile index stayed statistically unchanged up to 75% OCC (Table 2) and bulk/caliper remained consistent across this range (Table 4), the TSA results below were interpreted as reflecting contributions from both fiber/network-level characteristics and sheet structure.

Table 6 shows the softness of creped handsheets, measured with the Tissue Softness Analyzer. The surface softness parameter TS7, which captures high-frequency fiber vibration, decreased with more OCC and showed a stepwise response with a common B/C plateau. Sample A (100% BRF) had the highest TS7 at 24.3 dB, indicating greater surface roughness. Samples B (75% BRF/25% OCC) and C (50% BRF/50% OCC) had similar TS7 values of 18.5 and 18.4 dB, despite a 25% OCC difference. Samples D (25% BRF/75% OCC) and E (100% OCC) had lower values of 15.2 and 14.9 dB, respectively, with each forming its own statistically distinct group. Thus, the TS7 response was stepwise: Samples B and C formed a common intermediate group, while Sample D separated from the B/C group and Sample E separated further from Sample D. Tukey HSD confirmed four distinct TS7 groups (A; B/C; D; E). This TS7 gradient occurred mainly within the furnish range (Samples A to D), where bulk and caliper were statistically indistinguishable, indicating that the TS7 response changed despite limited variation in macroscopic bulk and caliper across these furnishes. Fiber-level morphology may contribute to the TS7 response, while bulk effects became more apparent at Sample E. This matches sensory results. In the sensory panel, Sample A ranked as the least-soft furnish, but pairwise significance in the full-panel analysis was limited to comparisons of A with Samples B and E, consistent with contributions from fiber/network-level factors in addition to bulk and caliper.

TS750, which is associated with the sheet’s low-frequency structural stiffness measured by the TSA, decreased steadily along the furnish gradient: 140.5 dB (Sample A), 99.7 dB (Sample B), 81.8 dB (Sample C), 62.2 dB (Sample D), and 50.2 dB (Sample E). This gradual decline across five levels contrasts with the two-tier pattern observed in bulk measurements, which remained statistically similar across Samples A to D before increasing sharply at Sample E (Table 4). The partial decoupling between TS750 and bulk indicates that TS750 is not merely a proxy for specific volume; rather, its sensitivity in the A to D range, where sheet structure remains stable, suggests it also detects fiber-level network stiffness and bonding changes associated with progressive OCC substitution. Only at Sample E, where both bulk and caliper experience a discrete structural transition, does the bulk response appear to dominate TS750, as the sharp rise in specific volume coincides with the sheet’s lowest resistance to low-frequency compression from the TSA blades. ANOVA with Tukey HSD confirmed three statistically distinct groups (A alone; B and C intermediate; D and E lowest; p < 0.001), offering finer resolution than the two-tier bulk pattern and indicating that TS750 responds to more than just architectural changes.

Table 6. TSA Instrumental Softness Parameters for Creped Handsheets (mean ± SD, Sled-Contact Surface)

TSA Full Acoustic Spectrum Analysis

Beyond the scalar TS7 and TS750 parameters, the full acoustic spectra (Fig. 3) recorded by the TSA provided additional insight into how OCC substitution altered the vibrational response of the creped sheet structure across the frequency domain. All five furnishes exhibited a characteristic unimodal spectral profile with a single dominant peak in the low-frequency range, corresponding to the compressional resonance of the sheet under blade loading. However, the spectrum differed markedly in peak amplitude, peak frequency, and spectral shape in ways that scalar parameters alone did not fully capture.

Peak amplitude decreased with increasing OCC content. Sample A (100% BRF) exhibited the highest low-frequency peak response of 139.8 dB at approximately 200 Hz, whereas Sample E (100% OCC) exhibited a substantially lower peak response of 49.9 dB at approximately 300 Hz. The progressive shift of the spectral peak from 200 Hz in Samples A and B to 250 Hz in Samples C and D, and to 300 Hz in Sample E, indicated that increasing OCC content shifted the dominant peak toward higher frequencies of the sheet under compressive blade loading, consistent with the concurrent increase in specific volume and the reduction in sheet stiffness associated with the higher-bulk, less-bonded OCC fiber network.

Fig. 3. Mean TSA acoustic spectra (0 to 8,000 Hz) for creped BRF/OCC handsheets (Samples A–E). Shaded bands represent ±1 standard deviation. Dashed vertical lines indicate the TS750 peak frequency (~200–250 Hz) and the TS7 readout frequency (6,500 Hz). Sled-contact surface, TP II algorithm; n = 5 per furnish.

The difference between the low-frequency peak amplitude and the TS7-region response at approximately 6,500 Hz decreased from 118.1 dB for Sample A to 37.7 dB for Sample E, indicating that the relative dominance of the low-frequency response decreased with increasing OCC content. This result suggests that the acoustic signature of BRF-rich sheets is more strongly dominated by the low-frequency structural response, whereas the relative contribution of the high-frequency surface response becomes greater as OCC content increases. Inspection of the full spectra showed that Sample C (50% BRF / 50% OCC) maintained a broader, more diffuse low-frequency response than the other furnishes. This broader response suggests that its creped sheet may exhibit greater heterogeneity in compressional compliance. This could result from the interaction between BRF and OCC fibers at equal blend ratios. This spectral anomaly aligns with Sample C’s unexpectedly low ranking in the sensory panel compared to Sample D. Importantly, this spectral heterogeneity occurs even though Sample C’s bulk and caliper were statistically similar to those of Samples A, B, and D (Table 4), suggesting that the anomaly reflects a unique compressional response not captured by bulk or caliper measurements.

Taken together, the spectral analyses revealed that OCC substitution produced a systematic depression and rightward shift of the dominant low-frequency resonance peak, accompanied by a relative increase in the contribution of high-frequency surface vibration. These spectral changes were consistent with the reduction in TS7 and TS750 values across the blend gradient but provided additional mechanistic details that scalar parameters alone did not convey.

The numerically integrated spectral-response (AUC) values for all five furnishes are presented in Table 7. The AUC decreased from 287,000 dB·Hz for Sample A to 233,000 dB·Hz (B), 199,000 dB·Hz (C), 158,000 dB·Hz (D), and 117,000 dB·Hz for Sample E, a total reduction of 59.1% across the blend gradient. The furnish means followed a strictly ordered sequence, providing a complementary descriptive measure of the overall TSA spectral response across the BRF-to-OCC substitution gradient. Because AUC was calculated from the mean spectrum for each furnish, these values were not used for replicate-level inferential comparisons among furnishes.

Table 7. Area under the Full TSA Acoustic Spectrum (AUC, 0–23,950 Hz) for Each Furnish, Computed by Numerical Integration of the Mean dB Curve

Sensory Panel Quality

Panel quality metrics indicated that the forced-choice sensory evaluation data were reliable and free from systematic bias. The mean panelist consistency, assessed by intransitive triad analysis, was 86.9%, indicating that most panelists made transitive and internally coherent choices across the 10 pairwise comparisons (full derivation in Supplementary Information, Section S7.1). Consistency values above 75 to 80% are generally considered acceptable for naïve panelists evaluating unlabeled tactile substrates (Siegel and Castellan 1988), and the mean value of 86.9% observed here indicates that the panel performed well above this threshold. A two-tailed binomial test applied to the aggregate right-hand win rate (53.5%; n = 260 comparisons) yielded p = 0.291, confirming that no significant positional hand bias was present at the α = 0.05 level (Supplementary Information, Section S7.4). Kendall’s coefficient of concordance for the full panel was W = 0.085 (χ²(4) = 8.86, p = 0.065; Supplementary Information, Section S7.3). Although this value did not reach conventional statistical significance and reflects low overall agreement in the complete ranking, this outcome is expected and does not indicate poor panel performance. Kendall’s W measures agreement across the full five-sample rank order simultaneously; when four of the five samples are perceptually close to one another, panelists will disagree on their relative ordering even while consistently identifying the outlier sample. In this context, the low W reflects the inherent difficulty of the discrimination task rather than panel unreliability and does not compromise the validity of the individual pairwise comparisons, which are evaluated independently and are not contingent on global rank agreement. Four panelists were excluded from the cleaned-panel analyses based on the outlier criteria described in the Data Analysis section (Supplementary Information, Section S7.2). The consistency and bias-check results collectively supported the use of Bradley–Terry modeling, Thurstone scaling, and pairwise binomial sign tests for characterizing the softness ranking among the five furnishes.

Sensory Panel: Softness Ranking

The raw win rates provided an initial indication of the perceived softness hierarchy. Samples E and B obtained the highest win rates (59.6% each), followed by Sample D (57.7%), Sample C (46.2%), and Sample A (26.9%). Sample A was the only furnish to fall substantially below the theoretical chance level of 50%, consistent with a systematic preference against the 100% BRF benchmark across the majority of pairwise comparisons in which it was evaluated.

Bradley–Terry modeling of the full-panel win matrix yielded latent softness strength scores of 0.083 for Sample A, 0.163 for Sample C, 0.241 for Sample D, 0.257 for Sample B, and 0.257 for Sample E, placing Samples B and E in a tied first position, with Sample D second, Sample C third, and Sample A last. After removal of the four panelists meeting the prespecified sensitivity-screening criterion, the scores were 0.069 (A), 0.141 (C), 0.226 (D), 0.254 (B), and 0.311 (E), with the rank order unchanged. Sample E’s score increased from 0.257 to 0.311 in this secondary sensitivity analysis; however, the reduced-panel results were not treated as the primary inferential result. The full win matrix, iterative update procedure, and convergence trace are provided in Supplementary Information, Section S3.

Thurstone Case V scaling confirmed the rank order and additionally quantified the perceptual distances between samples on a continuous interval scale. In the full panel analysis, scale values were −0.656 for Sample A, −0.101 for Sample C, +0.196 for Sample D, +0.261 for Sample B, and +0.300 for Sample E. In the cleaned-panel analysis, these values shifted to −0.810 (A), −0.185 (C), +0.176 (D), +0.281 (B), and +0.537 (E). In both analyses, Sample A had the lowest Thurstone scale value, with a separation of more than 0.55 scale units from the nearest sample (C), placing Sample A at the bottom of the model-derived softness scale. Within the OCC-containing group, Sample C consistently showed a negative scale value, indicating that its perceived softness was below the panel mean, whereas Samples B, D, and E all scored positively. The full z-score matrix, column means, and scale value derivation are provided in Supplementary Information, Section S4.

Pairwise binomial sign tests identified which differences were statistically significant at α = 0.05. For the full panel (n = 26), two pairs were significant: Sample E was preferred over Sample A with 88.5% agreement (p = 0.0001), and Sample B was preferred over Sample A with 76.9% agreement (p = 0.009). No significant differences were detected among any of the OCC-containing furnishes (Samples B, C, D, and E), with the closest non-significant result being B vs. C at 69.2% agreement (p = 0.076). For the cleaned panel (n = 22), the A vs. E comparison strengthened to 95.5% agreement (p < 0.0001), and a new significant difference emerged between Sample B and Sample C (77.3% agreement, p = 0.017), indicating that the 75% BRF / 25% OCC blend was perceived as softer than the 50% BRF / 50% OCC blend once the less consistent panelists were excluded. The A vs. B pair narrowly failed to reach significance in the cleaned panel (72.7% agreement, p = 0.053). The prevalence of non-significant results among pairs with win-rate differences in the 57 to 69% range was attributable to insufficient statistical power: post-hoc power calculations indicated that 51 panelists would have been required to achieve 80% power for a B vs. C type difference (69% agreement), and 83 to 147 panelists for the A vs. C and A vs. D comparisons (62 to 65% agreement). The current panel size of 26 was adequate only for large softness differences, such as A vs. E (power = 99.3% at n = 26). Complete p-values for all 10 pairs across both panels and full power calculations are provided in Supplementary Information, Sections S5 and S6.

Taken together, the win rates, Bradley–Terry scores, and Thurstone scale values showed a consistent directional ranking of E ≥ B > D > C > A. However, pairwise significance was limited: in the full-panel analysis, only A vs. B and A vs. E were significant, while the cleaned sensitivity analysis additionally identified B vs. C as significant. Most pairwise differences among the OCC-containing furnishes were not statistically significant at the current panel size; therefore, the sensory results are interpreted primarily as directional trends rather than definitive separation among all furnish levels.

Correlation between TSA and Sensory Panel

The degree of agreement between the TSA instrumental rankings and the sensory panel softness rankings is summarized in Table 8, which aligns the TS7 rank, TS750 rank, AUC rank, win rate, Bradley–Terry score, Thurstone scale value, and panel rank for each furnish alongside an assessment of concordance between the methods.

Table 8 shows directional agreement between the TSA metrics and the sensory panel at the endpoints: Sample E (100% OCC) ranked softest instrumentally, whereas Sample A (100% BRF) ranked least soft, and the sensory panel showed the same directional endpoint ordering. All three TSA metrics ranked the five furnishes in the order E (1^st) > D (2^nd) > C (3^rd) > B (4^th) > A (5^th), reflecting a systematic decrease in the integrated spectral response as OCC content increased. The sensory panel, however, ranked Sample B (75% BRF / 25% OCC) 2^nd rather than 4^th, placing it above Samples C and D; in the cleaned sensitivity analysis, B > C reached significance (p = 0.017). Because most pairwise differences among the OCC-containing furnishes were not significant in the full-panel analysis, the agreement is interpreted as a directional trend rather than evidence of complete sensory separation across all furnish levels.

Sample A (100% BRF) occupied the least-soft position in the TSA rankings and also ranked last in the sensory models, while Sample E showed agreement between the AUC and panel endpoint ranks. These results indicate directional consistency between instrumental and sensory responses at the blend endpoints, but they do not establish statistically significant sensory separation among all OCC-containing furnish levels.

Table 8 also highlights the different relative positions of Samples B, C, and D. The TSA ranked these three samples as D (2^nd) > C (3^rd) > B (4^th), reflecting a decrease in the integrated spectral-response index with increasing OCC content.

Table 8. Comparison of TSA Instrumental Rankings and Sensory Panel Softness Rankings across the Five BRF/OCC Furnishes

The sensory panel instead ranked Sample B 2^nd, above both C and D. This suggests that at 25% OCC substitution, the sensory response was not fully captured by the full-spectrum AUC ranking, potentially reflecting surface characteristics that are weighted differently by handfeel than by the integrated TSA response. The panel also ranked Sample C below Sample D, whereas the TSA metrics placed C above D, providing another indication that instrumental and sensory rankings were not identical across the intermediate blends.

To determine the individual effects of sheet architecture and fiber bonding on the TSA response, both bulk and tensile index were analyzed together as potential structural factors influencing TS7 and TS750. Bulk exhibited a stronger correlation with these parameters (TS750: r = -0.75; TS7: r = -0.67) compared to tensile index (TS750: r = 0.44; TS7: r = 0.36). A combined model incorporating both bulk and tensile index explained most of the variation in TS750 (R² = 0.97) and TS7 (R² = 0.88). This suggests that the TSA-measured softness is primarily controlled by bulk, with tensile index—and consequently, inter-fiber bonding—playing a secondary, complementary role, especially evident at the 100% OCC level where bulk increased sharply and tensile index decreased.

Correlation between Fiber Morphology, Paper Properties, and Softness Results

To analyze how furnish morphology relates to tissue performance, Pearson correlation coefficients were calculated between each of the seven FQA descriptors and the six tissue outcome variables across five different furnishes (Fig. 4).

Fig. 4. Pearson correlation coefficients (r) measure the relationship between fiber quality analyzer (FQA) morphological descriptors and tissue handsheet outcome variables across five BRF/OCC furnish blends (Samples A–E). Rows list FQA descriptors: length-weighted mean fiber length (Lw), ISO coarseness, fiber wall thickness (FWT), mean fiber width, curl index, kink index, and length-weighted fines content (LW Fines). Columns include tissue outcome variables: surface softness (TS7, dB), structural stiffness (TS750, dB), integrated TSA spectral-response index (AUC, dB·Hz), specific volume (Bulk, cm³/g), tensile index (kN·m/kg), and Thurstone perceived softness scale. Color intensity indicates the r value, ranging from −1 (red, strong negative correlation) to +1 (blue, strong positive correlation). Asterisks (*) mark correlations significant at α = 0.05 (two-tailed Pearson test, n = 5). BRF = bleached recycled fiber; OCC = old corrugated container.

The correlation matrix revealed a distinct macroscopic pattern: FQA descriptors exhibited strong, statistically significant correlations with instrumental acoustic softness (AUC and TS750) and specific volume. However, they did not show statistically significant correlations with tensile index or perceived softness as measured by Thurstone scaling.

Among all pairwise relationships, correlations with the integrated spectral response (AUC) were the strongest and most consistent. The kink index emerged as the strongest correlation with AUC (r = −0.987, R² = 0.975, p = 0.002), followed by Lw (r = −0.982, R² = 0.964, p = 0.003), coarseness (r = −0.977, R² = 0.954, p = 0.004), and width (r = −0.967, R² = 0.934, p = 0.007). Because the correlations are based on only five furnish-level means along a single BRF-to-OCC substitution gradient, they should be interpreted with caution and not as independent predictive or causal relationships. The strong co-variation among OCC content and fiber-morphology descriptors limits attribution to any single morphological parameter. This also suggests that AUC is a more morphologically sensitive softness descriptor than the scalar TS7 or TS750 readings alone. TS750 also showed similarly strong and significant correlations with kink index (r = 0.950, p = 0.013), Lw (r = −0.934, p = 0.020), coarseness (r = −0.928, p = 0.023), and width (r = −0.912, p = 0.031), indicating that bulk structural stiffness co-varied strongly with the morphological gradient across the five furnishes.

Correlations with TS7 typically aligned with those for TS750 and AUC, indicating that fibers that are longer, coarser, and straighter generally exhibit lower TS7 values. However, this pattern mainly results from the increased OCC in the furnish. Among the FQA descriptors, only the kink index showed a statistically significant correlation with TS7 (r = 0.905, p = 0.035). The weak and marginally significant links with TS7 align with its stepwise response with a B/C plateau, where Samples B and C had nearly identical TS7 values despite a 25-percentage-point difference in OCC content. This plateau effect restricts the range available for linear regression and somewhat separates the TS7 response from the overall morphological gradient.

Bulk was strongly and significantly correlated to all FQA descriptors related to fiber coarseness and cell wall dimensions, including FWT (r = 0.927, p = 0.024), Lw (r = 0.926, p = 0.024), coarseness (r = 0.915, p = 0.030), and width (r = 0.921, p = 0.027). Additionally, curl index and kink index were significantly correlated with bulk (r = −0.904, p = 0.035 and r = −0.889, p = 0.043, respectively).

No FQA descriptor showed a statistically significant correlation with the tensile index (best: FWT, r = −0.747, p = 0.147) or Thurstone perceived softness scale (best: kink index, r = −0.769, p = 0.128). This absence of significant FQA–tensile relationships matches the tensile index pattern shown in Table 2, where Samples B, C, and D had similar values that did not correlate with OCC content. Although FQA–TSA relationships were strong, the weak FQA–Thurstone correlations imply that the panel’s ranking (E > B > D > C > A) was not fully aligned with individual fiber-morphology descriptors. Perceived softness also depends on sheet-level surface effects not captured by FQA.

Screening Techno-economic and Monte Carlo Results

The cost distributions for both routes showed considerable overlap: treated OCC had a median headbox-ready cost of USD 517/odt (P10–P90: USD 425–621), while conventional high-grade recovered paper was USD 598/odt (P10–P90: USD 510–696) (Table 9; Fig. 5a). The median economic benefit, derived from subtracting treated OCC costs from conventional ones, was USD 80/odt, with a P10 to P90 range of USD -54 to 214/odt. Treated OCC was less costly in 77.7% of 100,000 trials (Fig. 5b). This model’s results depend on the distributions and the assumption of independence rather than actual outcomes. The significant overlap suggests that feedstock price alone does not determine the preferred route; other factors such as yield, treatment costs, processing difficulty, reject generation, and fixed costs also play crucial roles in the comparison.

Table 9. Summary Statistics from a Monte Carlo Simulation (100,000 Trials) Comparing Treated-OCC and Conventional RCP Routes

Fig. 5. Monte Carlo simulation (n = 100,000) compares the techno-economic aspects of treated OCC and high-grade recovered paper (RCP). (a) Fiber cost distributions for treated-OCC (blue solid) and RCP (orange dashed). (b) Economic benefit of treated OCC, calculated as conventional cost minus treated-OCC cost; the vertical dashed line shows break-even, and the dotted line shows median advantage. Treated OCC was cheaper in 77.7% of trials.

The median break-even OCC price was USD 252/odt (P10–P90: USD 167–348), representing the highest price at which treated OCC aligns with conventional-route costs under the sampled conditions. This threshold serves as a practical decision point rather than an unconditional cost claim: below it, treated OCC is likely to be more cost-effective; above it, factors like yield loss or processing expenses may negate the benefit of using cheaper feedstock. In this refined model, no trials indicated a negative break-even price, suggesting that treated OCC remained competitive at some positive purchase price in all 100,000 simulations. However, the wide range means the value cannot be generalized across mills without specific data on yield, chemicals, energy, and disposal costs.

Sensitivity analysis indicates that the main factors affecting route competitiveness are the delivered prices and yields for both routes. This is followed by the fixed processing allowance of the OCC route and the conversion and fixed costs associated with the conventional route (see Fig. 6). Yield has a disproportionately large, nonlinear impact because lower yields increase both input needs and disposal volumes simultaneously. This makes yield more important than an equivalent uncertainty in individual chemical prices or reject-disposal costs. These results depend on the input ranges used and assume independent inputs; in reality, cheaper recovered paper might contain more contaminants (lower yield), and more intensive treatment could increase both costs and yield recovery. Therefore, this analysis should be viewed as highlighting data-collection priorities rather than establishing universal economic drivers.

Fig. 6. Global sensitivity of the modeled economic advantage of treated OCC to six key uncertain inputs, shown as Spearman rank correlations. Positive (blue) means higher input increases advantage; negative (red) means the opposite. Coefficients depend on input ranges and independence assumptions, not causal effects.

Integration of Technical Performance and Furnish-level Economics

Economic outcomes should be interpreted alongside the observed trade-offs between softness and strength from experiments (Table 3; Fig. 2). The tensile index stayed the same or was slightly higher at 25 to 75% OCC substitution (7.77, 7.99, and 7.58 kN·m/kg for 25%, 50%, and 75% OCC, respectively) compared to the 100% BRF control (7.21 kN·m/kg). Meanwhile, TS7 and TS750 steadily decreased across this range (Table 6). At 100% OCC substitution, bulk reached its peak while TS7 and TS750 hit their lowest points; however, the tensile index sharply dropped to 4.91 kN·m/kg, a 31.9% reduction from the control (Table 2). Neither the lowest furnish cost nor the highest OCC fraction alone indicated an optimum: furnish cost decreases steadily from USD 598/odt at 0% OCC to USD 517/odt at 100% OCC (Table 9), while tensile index only significantly declined between 75% and 100% substitution. Consequently, the cost and strength curves only diverge notably at the highest substitution level (Fig. 7).

An economically optimal OCC fraction requires clear constraints on factors such as tensile index, stretch, TEA, bulk, absorbency, brightness, drainage, drying demand, sheet-break frequency, machine speed, converting yield, and saleable output. Because the blend-cost relationship is roughly linear and does not account for nonlinear production penalties, it cannot determine a jointly optimal OCC percentage. As a result, the 25 to 75% range remains a viable candidate, being the range where tensile index stayed at or above the control level while furnish cost decreased by USD 22 to 62 per odt (Table 9).

Fig. 7. The median furnish-level cost, derived from Monte Carlo simulation with error bars indicating the P10 to P90 range, and the measured tensile index across five BRF/OCC blend ratios (Samples A–E) are shown. Furnish cost decreases steadily as OCC substitution increases, while the tensile index stays statistically similar to the 100% BRF control (see Table 2) up to 75% OCC, then drops by 31.9% at 100% OCC.

The median USD 80/odt modeled cost advantage represents furnish-level economic headroom that could be reduced or eliminated by additional costs and performance penalties not considered in the present screening-level model. The notable drop in tensile index at 100% OCC substitution is an example of such an excluded penalty: since the screening model treats yield as a mass-recovery factor independent of sheet strength, the apparent cost benefit at 100% OCC may overestimate actual value if strength-related machine or converting losses occur. Conversely, benefits not included in the model, such as reduced dependence on premium recovered-paper grades or improved utilization of local OCC, could strengthen the case at moderate substitution levels. These effects are mill-specific and were not quantified here.

CONCLUSIONS

1. Replacing bleached recycled fiber (BRF) with bleached old corrugated container (OCC) in creped tissue increased bulk while maintaining a similar tensile index up to 75% OCC substitution.
2. At 100% OCC substitution, tensile index decreased by 31.9%, consistent with accompanying increases in bulk and caliper, indicating a strength penalty at full replacement.
3. OCC substitution was associated with changes in mechanical properties and TSA metrics, while the sensory panel showed the same directional endpoint trend; both approaches identified 100% BRF as the least soft furnish.
4. Across the five furnishes, TSA response showed strong associations with several fiber-morphology descriptors; however, these relationships should be interpreted with caution because the variables co-varied along the same BRF-to-OCC substitution gradient.
5. Monte Carlo analysis showed that treated OCC was lower cost than conventional recovered paper in most simulated market conditions, with a median modeled furnish-level cost advantage of USD 80/odt and a break‑even price of approximately USD 252/odt; these advantages depended primarily on yield and delivered price rather than feedstock cost alone.
6. The most favorable balance between modeled furnish-level cost advantage and tensile performance occurred at 25 to 75% OCC substitution, corresponding to USD 22 to 62/odt in modeled furnish-level cost advantage without detectable tensile-index loss, whereas at 100% OCC the trade-off between modeled cost advantage and tensile performance indicates that the preferred furnish level will depend on mill-specific strength constraints and converting requirements.

ACKNOWLEDGMENTS

The authors are grateful for the support of the Sustainable and Alternative Fibers Initiative (SAFI) members.

Conflict of Interest

The authors declare no conflict of interest

Use of Generative AI

The authors used an AI‑based assistant (Perplexity) for language editing and grammar correction during manuscript preparation; all scientific content, data analysis, and interpretation were conceived and validated by the authors.

REFERENCES CITED

Abdullah, R., Ishak, C., Kadir, W., and Bakar, R. (2015). “Characterization and feasibility assessment of recycled paper mill sludges for land application in relation to the environment,” Int. J. Environ. Res. Public Health 12(8), 9314-9329. https://doi.org/10.3390/ijerph120809314

Berg, P., and Lingqvist, O. (2019). “Pulp, paper, and packaging in the next decade: Transformational change,” (https://www.mckinsey.com/industries/packaging-and-paper/our-insights/pulp-paper-and-packaging-in-the-next-decade-transformational-change#/), accessed June 28, 2026.

Bradley, R. A., and Terry, M. E. (1952). “Rank analysis of incomplete block designs: I. The method of paired comparisons,” Biometrika 39(3-4), 324-345. https://doi.org/10.1093/biomet/39.3-4.324

Christopher, M. (2020). “Focus on: Optimising deink plant quality and yield,” (https://www.tissueworldmagazine.com/special-features/focus-on/focus-on-optimising-deink-plant-quality-and-yield/), accessed August 20, 2026.

D’Olivo, I. (2022). “5 essential properties of a soft and bulky tissue paper,” (https://www.compendiumpaperasia.com/index.php/cpa/article/view/9), accessed June 28, 2026.

de Assis, T., Pawlak, J., Pal, L., Jameel, H., Venditti, R., Reisinger, L. W., Kavalew, D., and Gonzalez, R. W. (2019). “Comparison of wood and non-wood market pulps for tissue paper application,” BioResources 14(3), 6781-6810. https://doi.org/10.15376/biores.14.3.6781-6810

de Assis, T., Reisinger, L. W., Pal, L., Pawlak, J., Jameel, H., and Gonzalez, R. W. (2018). “Understanding the effect of machine technology and cellulosic fibers on tissue properties – A review,” BioResources 13(2), 4593-4629. https://doi.org/10.15376/biores.13.2.4593-4629

EMTEC Electronic GmbH (2025). “Technical basics of tissue softness measurement,” (https://industrialphysics.com/knowledgebase/articles/technical-basics-of-tissue-softness-measurement/), accessed August 20, 2026.

Hämäläinen, P., Hallbäck, N., and Barbier, C. (2016). “Development and evaluation of a high-speed creping simulator for tissue,” Nord. Pulp Paper Res. J. 31(3), 448-458. https://doi.org/10.3183/npprj-2016-31-03-p448-458

Hollmark, H., and Ampulski, R. S. (2004). “Measurement of tissue paper softness: A literature review,” Nord. Pulp Paper Res. J. 19(3), 345-353. https://doi.org/10.3183/npprj-2004-19-03-p345-353

ISO 12625-3 (2014). “Tissue paper and tissue products – Part 3: Determination of thickness, bulking thickness and apparent bulk density,” International Organization for Standardization, Geneva, Switzerland.

Kendall, M. G. (1975). Rank Correlation Methods, 4^th Ed., Griffin, London, UK.

Moore, B., and Cornish, S. (2020). “Recovered paper supply – already in decline, hit hard by the pandemic,” (https://www.tissueworldmagazine.com/departments/exitissues/recovered-paper-supply-already-in-decline-hit-hard-by-the-pandemic/), accessed June 28, 2026.

Morris, H. (2018). “Analysis: How supply and costs are changing across the world’s recovered paper markets?” (https://www.tissueworldmagazine.com/departments/exitissues/analysis-how-supply-and-costs-are-changing-across-the-worlds-recovered-paper-markets/), accessed August 2, 2026.

Padley, I. (2022). “The basics of creping in the tissuemaking process,” (https://www.tissuestory.com/2022/04/26/the-basics-of-creping/), accessed June 28, 2026.

Pawlak, J. J., Frazier, R., Vera, R. E., Wang, Y., and Gonzalez, R. (2022). “Review: The softness of hygiene tissue,” BioResources 17(2), 3509-3550. https://doi.org/10.15376/biores.17.2.pawlak

Prinz, M., Zollner-Croll, H., Bauer, W., and Meltzer, F. P. (2021). “A new method to evaluate the hand-feel potential of furnishes for tissue paper products,” BioResources 16(1), 684-697. https://doi.org/10.15376/biores.16.1.684-697

Quinn, M. (2024). “MRF and mill advances expected to fuel North American recovered paper markets,” (https://www.wastedive.com/news/recovered-paper-outlook-occ-mixed-paper-fastmarkets/726956/), accessed June 28, 2026.

Reynolds, A. M., Franco-Ocampo, J., Hubbe, M. A., Park, S., and Venditti, R. A. (2025). “Fractionation by classification as a prerequisite to the purification of old corrugated container material for preparing dissolving grade pulp,” Cellulose 32, 6209-6228. https://doi.org/10.1007/s10570-025-06530-x

Siegel, S., and Castellan, N. J. (1988). Nonparametric Statistics for the Behavioral Sciences, 2^nd Ed., McGraw-Hill, New York, NY, USA.

Tanoh, R., Nikiema, J., Asiedu, Z., Jayathilake, N., and Cofie, O. (2022). “The contribution of tipping fees to the operation, maintenance, and management of fecal sludge treatment plants: The case of Ghana,” J. Environ. Manage. 303, article 114125. https://doi.org/10.1016/j.jenvman.2021.114125

TAPPI T 410 om-23 (2023). “Grammage of paper and paperboard (weight per unit area),” TAPPI Press, Atlanta, GA, USA.

TAPPI T 581 om-22 (2022). “Dry tensile properties of paper towel and tissue products (using constant rate of elongation apparatus),” TAPPI Press, Atlanta, GA, USA.

TAPPI T 402 sp-21 (2021). “Standard conditioning and testing atmospheres for paper, board, pulp handsheets, and related products,” TAPPI Press, Atlanta, GA, USA.

Thurstone, L. L. (1927). “A law of comparative judgment,” Psychol. Rev. 34, 273-286. https://doi.org/10.1037/h0070288

TissueOnline (2023). “Tissue sector faces shortage of recyclable paper supply,” (https://tissueonlinenorthamerica.com/tissue-sector-faces-shortage-of-recyclable-paper-supply/), accessed June 28, 2026.

Urdaneta, F., Kumar, R., Marquez, R., Vera, R. E., Franco, J., Urdaneta, I., Saloni, D., Venditti, R. A., Pawlak, J. J., Jameel, H., and Gonzalez, R. W. (2024). “Evaluating chemi-mechanical pulping processes of agricultural residues: High-yield pulps from wheat straw for fiber-based bioproducts,” Ind. Crops Prod. 221, article 119379. https://doi.org/10.1016/j.indcrop.2024.119379

Urdaneta, F., Meza, L., Marquez, R., Franco, J., Vera, R. E., Urdaneta, I., Parrish, C., Saloni, D., Pawlak, J. J., Jameel, H., and Gonzalez, R. W. (2025). “Comparing the impacts of wood and APMP non-wood fibers on the properties of hygiene tissue paper,” Cellulose 32, 5021-5034. https://doi.org/10.1007/s10570-025-06529-4

Urdaneta, I., Vera, R. E., Marquez, R., Vivas, K. A., Urdaneta, F., Gongora, S., Frazier, R., Franco, J., Cordoba, Y., Jameel, H., Viitala, R., Harbec, G., Carette, J., Pawlak, J. J., and Gonzalez, R. W. (2025). “Assessing and predicting the softness of hygiene tissue containing alternative fibers,” Cellulose 32(2), 1171-1199. https://doi.org/10.1007/s10570-024-06289-7

U.S. Environmental Protection Agency (2000). “Mills that produce tissue products containing recovered paper,” Washington, DC, USA.

Veijalainen, S. (2022). “How fiber and energy mixes impact containerboard price,” (https://www.fastmarkets.com/insights/how-fiber-and-energy-mixes-impact-containerboard-price/), accessed August 2, 2026.

Vieira, J. C., Fiadeiro, P. T., and Costa, A. P. (2023). “Converting operations impact on tissue paper product properties – A review,” BioResources 18(1), 2303-2326. https://doi.org/10.15376/biores.18.1.Vieira

Wang, Y., de Assis, T., Zambrano, F., Pal, L., Venditti, R., Dasmohapatra, S., Pawlak, J., and Gonzalez, R. (2019). “Relationship between human perception of softness and instrument measurements,” BioResources 14(1), 780-795. https://doi.org/10.15376/biores.14.1.780-795

Workman, M. (2024). “US OCC, DLK, mixed paper pricing runs show no signs of slowdown,” (https://www.fastmarkets.com/insights/us-occ-dlk-mixed-paper-pricing-runs-show-no-signs-of-slowdown/), accessed August 2, 2026.

World Paper Mill (2026). “US printing-writing paper shipments decrease November.”

Zambrano, F., Marquez, R., Jameel, H., Venditti, R. A., and Gonzalez, R. W. (2021). “Upcycling strategies for old corrugated containerboard to attain high-performance tissue paper: A viable answer to the packaging waste generation dilemma,” Resour. Conserv. Recycl. 175, article 105854. https://doi.org/10.1016/j.resconrec.2021.105854

Zambrano, F., Marquez, R., Vera, R. E., Jameel, H., Venditti, R. A., and Gonzalez, R. W. (2022). “Developing alternative, high-absorbency brown fibers: Tissue paper from upcycled corrugated packaging waste to meet new consumer trends,” ACS Sustain. Chem. Eng. 10, 13343-13356. https://doi.org/10.1021/acssuschemeng.2c03280

Article submitted: August 10, 2026; Peer review completed: September 12, 2026; Revised version received: September 12, 2026; Accepted: September 13, 2026; Published: September 28, 2026.

DOI: 10.15376/biores.21.4.11129-11166

APPENDIX

S1. AREA UNDER THE TSA ACOUSTIC SPECTRUM (AUC)

S1.1 Definition and Formula

The area under the full TSA response spectrum (AUC) was computed by trapezoidal numerical integration of the mean TSA response A(f) (in dB) over the complete frequency range recorded by the instrument (0 to 23,950 Hz). For each furnish, replicate spectra were first averaged point-by-point to obtain a mean spectrum, which was then integrated. Accordingly, the AUC is treated as a descriptive furnish-level spectral-response index rather than as a replicate-level inferential response variable. Because the TSA output is reported in dB, the resulting dB·Hz quantity should not be interpreted as physical acoustic energy. The trapezoidal rule approximates the integral as the sum of trapezoids formed between adjacent frequency points:

AUC = ∫₀²³’⁹⁵⁰ A(f) df ≈ Σᵢ ½ [A(fᵢ) + A(fᵢ₊₁)] × (fᵢ₊₁ − fᵢ) (S1)

where A(fᵢ) is the mean TSA response in dB at frequency fᵢ, and the summation runs over all adjacent frequency intervals recorded from 0 to 23,950 Hz. The AUC is expressed in units of dB·Hz and is used here only as an integrated spectral-response index. Calculations were implemented using numpy.trapezoid in Python 3.10.

S1.2 Worked Example — Sample A (100% BRF)

The first five and final trapezoidal steps for Sample A are shown below. Mean dB values are taken from the averaged replicate spectrum.

Table S1. Selected Trapezoidal Steps for Sample A AUC Calculation

Table S2. Full-spectrum AUC for All Five Furnishes

S2. WIN RATE

S2.1 Experimental Design

Softness was evaluated using a complete round-robin pairwise comparison design. With five furnishes (A–E), all C(5,2) = 10 possible pairs were evaluated: A–B, A–C, A–D, A–E, B–C, B–D, B–E, C–D, C–E, and D–E. Each pair was presented to all 26 panelists (full panel), yielding 26 × 10 = 260 total pairwise judgments. Because each sample appears in exactly four of the ten pairs, each sample received 4 × 26 = 104 individual judgments in total.

S2.2 Definition and Formula

The win rate Wᵢ for sample i is the percentage of the 104 total judgments in which that sample was chosen as the softer of the two:

Wᵢ = (votes won by sample i) / 104 × 100% (S2)

S2.3 Worked Example — Sample A

Sample A was involved in four pairs: A vs. B, A vs. C, A vs. D, and A vs. E. Across all 104 panelist judgments involving Sample A, panelists chose A as softer in 28 instances. Therefore:

W_A = 28 / 104 × 100% = 26.9% (S3)

Table S3. Win Rates for All Five Furnishes (full panel, n = 26)

S3. BRADLEY–TERRY PAIRED-COMPARISON MODEL

S3.1 Win Matrix

The Bradley–Terry model is fitted to the win matrix W, where entry W[i][j] is the number of panelists who preferred sample i over sample j. The win matrix for the full panel (n = 26) is shown in Table S4. The matrix is asymmetric: W[i][j] + W[j][i] = 26 for every off-diagonal pair (each panelist voted once per pair). Diagonal entries are zero by definition.

Table S4. Win Matrix W[i][j] = Number of Panelists Who Preferred Row Sample Over Column Sample (full panel, n = 26)

Row sum equals votes won (Table S3). Column sum equals votes lost. W[i][j] + W[j][i] = 26 for all pairs.

S3.2 Model Definition and Iterative Update Rule

The Bradley–Terry model (Bradley and Terry 1952) assigns a latent softness strength parameter πᵢ > 0 to each sample i, such that the probability of sample i being judged softer than sample j is:

P(i ≻ j) = πᵢ / (πᵢ + πⱼ) (S4)

Parameters are estimated by maximum likelihood through iterative scaling. At each iteration, the parameter for sample i is updated as:

πᵢ⁽ᵗ₊¹⁾ = Wᵢ / Σⱼ≠ᵢ [(Wᵢⱼ + Wⱼᵢ) / (πᵢ⁽ᵗ⁾ + πⱼ⁽ᵗ⁾)] (S5)

where Wᵢ = Σⱼ W[i][j] is the total number of wins for sample i, and the denominator sums the expected number of comparisons involving i weighted by the current parameter estimates. After each iteration, all parameters are renormalised to sum to 1. Convergence is declared when the maximum absolute change across all parameters falls below 10⁻⁶. For this dataset convergence was achieved in 21 iterations.

S3.3 Convergence Trace

Table S5. Bradley–Terry Iterative Scaling — Selected Iterations (full panel)

S3.4 Worked Example — Predictive Probabilities

Using the converged scores, the model-predicted probabilities for key comparisons are:

P(E ≻ A) = 0.257 / (0.257 + 0.082) = 0.758

P(B ≻ A) = 0.257 / (0.257 + 0.082) = 0.758

P(B ≻ C) = 0.257 / (0.257 + 0.163) = 0.612

Table S6. Bradley–Terry Latent Strength Scores — Full and Cleaned Panels

S4. THURSTONE CASE V SCALING

S4.1 Method

Thurstone Case V scaling (Thurstone 1927) converts the pairwise preference proportions into a continuous psychological interval scale. Let pᵢⱼ be the proportion of panelists who preferred sample i over sample j. Under Case V assumptions (equal and independent discriminal dispersions), each entry in the z-score matrix is:

zᵢⱼ = Φ⁻¹(pᵢⱼ) (S6)

Diagonal entries are set to 0. Proportions of exactly 0.500 yield z = 0. The raw scale value for sample j is the mean of column j, computed over the n−1 = 4 off-diagonal entries (the diagonal zero is excluded to avoid diluting the signal):

Scaleⱼ (raw) = (1/(n−1)) × Σᵢ≠ⱼ zᵢⱼ (S7)

Because a high column mean for sample j indicates that other samples are frequently preferred over j (i.e., j is perceived as less soft), the scale is negated so that higher values correspond to greater perceived softness. Finally, scale values are centered by subtracting their grand mean so that they sum to zero.

S4.2 Observed Proportions Matrix (Full Panel, n = 26)

Entry pᵢⱼ is the proportion of the 26 panelists who preferred sample i over sample j. Note that pᵢⱼ + pⱼᵢ = 1 for all i ≠ j.

Table S7. Pairwise Preference Proportions p(row preferred over column) — Full Panel

S4.3 z-Score Matrix

Each off-diagonal proportion is transformed to a z-score using the inverse normal (probit) function. For example: p(E ≻ A) = 0.885 → z = Φ⁻¹(0.885) = +1.200; p(A ≻ E) = 0.115 → z = Φ⁻¹(0.115) = −1.200. Note that zᵢⱼ = −zⱼᵢ, confirming internal consistency.

Table S8. Z-score Matrix Z (diagonal = 0; off-diagonal zᵢⱼ = Φ⁻¹(pᵢⱼ))

S4.4 Scale Value Derivation — Worked Example for Sample A

The column mean for Sample A (excluding diagonal) is the mean of the four off-diagonal entries in column A: (+0.736 + 0.292 + 0.396 + 1.200) / 4 = 2.624 / 4 = +0.656. After negation: −0.656. The grand mean of all five negated column means is (−0.656 + 0.261 − 0.101 + 0.196 + 0.300) / 5 = 0.000. Since the grand mean is already zero, no further centering adjustment is needed for this dataset. The final Thurstone scale values are therefore equal to the negated column means.

Table S9. Thurstone Case V Scale Values — Full and Cleaned Panels

S5. PAIRWISE BINOMIAL SIGN TESTS

S5.1 Test Procedure

For each of the 10 pairwise combinations, a two-tailed binomial sign test was applied. Under H₀, each panelist is equally likely to prefer either sample (p₀ = 0.50). Let k be the number of panelists who preferred the majority sample out of n total panelists. The two-tailed p-value is:

p = 2 × P(X ≥ k | n, 0.5) = 2 × Σᵥᵏ₌ᵏ C(n, m) × 0.5ⁿ (S8)

where X ~ Binomial(n, 0.5). For n = 26 the two-tailed rejection region at α = 0.05 is R = {k ≤ 7 or k ≥ 19}. For the cleaned panel (n = 22) the rejection region is R = {k ≤ 5 or k ≥ 17}.

S5.2 Worked Example — A vs. E (Full Panel, n = 26)

Sample E was preferred by 23 of 26 panelists (k = 23, 88.5%). Since k = 23 ≥ 19 it falls in the rejection region:

p = 2 × P(X ≥ 23 | 26, 0.5) = 2 × [C(26,23) + C(26,24)

+ C(26,25) + C(26,26)] × 0.5²⁶ = 0.0001 (S9)

Table S10. Pairwise Binomial Sign Tests — Full Panel (n = 26) and Cleaned Panel (n = 22)

S6. POST-HOC POWER ANALYSIS

S6.1 Method

Post-hoc power was calculated as the probability of rejecting H₀ given the observed agreement proportion p_obs and panel size n. Power equals the sum of binomial probabilities over the rejection region R = {k : p-value ≤ 0.05}:

Power = Σᵏ∈R P(X = k | n, p_obs), where X ~ Binomial(n, p_obs) (S10)

For n = 26, R = {0,1,2,3,4,5,6,7,19,20,21,22,23,24,25,26}. For n = 22, R = {0,1,2,3,4,5,17,18,19,20,21,22}. Sample sizes for 80% and 90% power were found by incrementing n until the calculated power exceeded each threshold.

S6.2 Worked Example — A vs. E (n = 26, p_obs = 0.885)

Using X ~ Binomial(26, 0.885) and R = {0–7, 19–26}, power is the sum of P(X = k | 26, 0.885) over all k ∈ R. The dominant contribution comes from the upper tail: P(X ≥ 19 | 26, 0.885) ≈ 0.993. The lower tail contribution P(X ≤ 7 | 26, 0.885) ≈ 0.000. Total power = 0.993.

Table S11. Post-hoc Power and Sample Size Requirements

S7. PANELIST QUALITY METRICS

S7.1 Intransitive Triad Analysis

With five samples evaluated via all 10 pairwise comparisons, the total number of three-sample triads per panelist is C(5,3) = 10. A triad {i, j, k} is intransitive (circular) if the panelist’s three choices form a cycle: i ≻ j, j ≻ k, and k ≻ i simultaneously. Such cycles indicate internal inconsistency. Individual consistency was calculated as:

Consistency (%) = [(T − d) / T] × 100 (S11)

where T = 10 is the total number of triads and d is the number of circular triads for that panelist. A panelist with d = 0 achieves 100% consistency; d = 4 gives 60% consistency. The mean across 26 panelists was 86.9% (mean d ≈ 1.3 circular triads per panelist). Consistency above 75–80% is considered acceptable for naïve panelists on unlabeled tactile substrates (Siegel and Castellan 1988).

S7.2 Outlier Identification Criterion

Each panelist’s individual ranking was compared with the aggregate ordering derived from the full-panel analyses. Kendall’s τ was used to quantify concordance between the panelist-level ordering and the aggregate ordering across the 10 possible sample pairs. The full 26-panelist dataset remained the primary sensory analysis; this screening was used only to construct a secondary sensitivity analysis.

τ = (n_concordant − n_discordant) / C(n, 2) (S12)

where n_concordant is the number of sample pairs ranked in the same order as the aggregate ordering, n_discordant is the number ranked in the opposite order, and C(n, 2) = 10 is the total number of pairwise comparisons. Kendall’s τ ≤ −0.2 was used as a pragmatic screening threshold because it represents more discordant than concordant pairwise rankings; with 10 comparisons, τ = −0.2 corresponds to approximately four concordant and six discordant comparisons. Four panelists met this criterion, with τ values from −0.2 to −0.4. This threshold was not treated as an independently validated sensory-outlier standard. Because the criterion was defined relative to the observed aggregate ordering, the reduced-panel results were interpreted only as a sensitivity analysis and not as the primary inferential result.

S7.3 Kendall’s Coefficient of Concordance (W)

Kendall’s W quantifies agreement among all k panelists in their complete rank ordering of n samples. Ranks are derived from each panelist’s individual win counts across the 10 pairwise comparisons: the sample winning the most comparisons for that panelist is ranked 1^st (softest), the sample with the second most wins is ranked 2^nd, and so on. Ties in win count are broken by the specific pair outcomes. The rank sum Rⱼ for sample j is the sum of that sample’s ranks across all k = 26 panelists. The expected rank sum under complete disagreement is:

Ā = k(n + 1) / 2 = 26 × 6 / 2 = 78

The deviation sum S is:

S = Σⱼ (Rⱼ − Ā)² (S13)

and Kendall’s W is:

W = 12S / [k²(n³ − n)] (S14)

W ranges from 0 (no agreement) to 1 (perfect agreement). The chi-squared test for H₀: W = 0 uses:

χ² = k(n − 1)W, df = n − 1 (S15)

For the full panel: k = 26, n = 5, W = 0.085. Therefore χ² = 26 × 4 × 0.085 = 8.84 (p = 0.065, df = 4). The non-significant result does not invalidate the panel; low W is expected and appropriate when most samples are perceptually similar and discrimination is difficult.

S7.4 Hand Position Bias Test

Each pairwise comparison was presented with one sample in the panelist’s left hand and one in the right hand. The aggregate right-hand win rate was 53.5% (139 right-hand wins out of 260 total comparisons). A two-tailed binomial test against H₀: p = 0.50 yielded p = 0.291. The 95% confidence interval for the true right-hand win probability is [47.3%, 59.7%], which encompasses 0.50, confirming no significant positional bias.

S8. References

Bradley, R. A., and Terry, M. E. (1952). “Rank analysis of incomplete block designs: I. The method of paired comparisons,” Biometrika 39(3-4), 324-345. https://doi.org/10.2307/2334029

Kendall, M. G. (1948). Rank Correlation Methods, Griffin, London, UK.

Siegel, S., and Castellan, N. J. (1988). Nonparametric Statistics for the Behavioral Sciences, 2^nd Ed., McGraw-Hill, New York, NY, USA.

Thurstone, L. L. (1927). “A law of comparative judgment,” Psychological Review 34(4), 273-286. https://doi.org/10.1037/h0070288