Modelling plasticiser loss inside closed environments

Plasticiser loss, a fundamental mechanism in the degradation of polymeric materials, leads to material instability and contamination of the environment. The process depends on environmental conditions, but the size of the container in which an object is housed and the thickness of the material also play key roles in the rate of plasticiser loss and the time at which equilibrium is reached. Understanding these dependencies provides valuable insight into the degradation of plastic museum artefacts inside enclosures as typified by museum storage and, more broadly, the deterioration of polymeric materials in closed environments. Migration of low molecular weight plasticisers, like diethyl phthalate or dimethyl phthalate, from plastics has been widely studied in accelerated ageing experiments at elevated temperatures and different airflows. Here, to investigate these effects, we modelled plasticiser loss in a stagnant environment inside an enclosure at room temperature. Our model is one-dimensional and describes loss through a two-phase transient diffusion process, between the solid plastic and the air. The comparison of numerical simulations to FTIR and 1 H NMR spectroscopic data from cellulose acetate samples plasticised with diethyl phthalate aged at T = 70 ◦ C and 50% relative humidity indicates that the model is appropriate for thin enclosed plastics. We applied the model over a range of diffusion coefficients [10 (cid:0) 21 – 10 (cid:0) 14 m 2 s (cid:0) 1 ] and partition coefficients [10 3 – 10 7 ] to represent different polymeric materials. Under the investigated scenarios, thin plastics tolerate a maximum total plasticiser loss of 10% and the timescales vary between 1 and 10 9 years. The model provides insight into the relationship between plasticiser loss, enclosure dimensions and material thickness for different plastics and can suggest how to improve packaging dimensions of thin plastic products to minimise loss, as well as future conservation strategies and guidelines for plastic museum artefacts.


Introduction
Plastics have become an integral part of modern life since their invention in the late 19th century, and their stability has been much discussed. Plastic stability is important in different fields to ensure safe operation, including industrial applications [1][2][3][4] and medical applications [5][6][7][8]. Furthermore, the long-term stability of plastics in a heritage context is significant, since nowadays plastic artefacts constitute a major part of museum collections [9].
Plasticisers are additives used during plastic manufacturing to decrease the glass transition temperature of the material and therefore enhance its flexibility [10]. Low molecular weight plasticisers like diethyl phthalate (DEP), dimethyl phthalate (DMP) or dibutyl phthalate (DBP) were historically used with cellulose acetate (CA) [11] in films [12][13][14][15] and artworks [16]. For example, CA Walt Disney animation cels have been found to incorporate all three plasticisers [15]. VOC analysis with SPME GC-MS fibres has also revealed emissions of both DEP and DMP from the CA-based Construction in Space 'Two Cones' by Naum Gabo in Tate's collection [16].
Additionally, plasticiser migration from plasticised polymers may change the surface texture with the deposition of crystals, bubbles [13] or sticky residues at the plastic surface, and, in terms of mechanical damage, leads to material warping, crazing, and ultimately cracking and embrittlement (Fig. 1). Regarding plastic artefacts in museums, these physical and chemical changes affect the stability but also the aesthetics of artefacts irreversibly, inhibiting their handling, displaying or application of treatments due to material vulnerability.
Degradation due to plasticiser loss depends not only on environmental factors, like temperature [33], relative humidity [34] and ventilation [38], but also on storage type, (its size and also whether it is sealed or open), storage material and material properties, such as plasticiser volatility and plastic thickness. However, there has been little investigation on the impact of the container size and plastic thickness in the rate of plasticiser loss and on the time that this process requires to reach equilibrium. These questions will be addressed in this paper.
Plasticiser loss has previously been investigated by using mathematical modelling based on diffusion and evaporation mechanisms [2,39,40] or convective mass transfer mechanism [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37] and by conducting accelerated ageing experiments at high temperatures and different airflow conditions [38,[41][42][43][44][45]. In this work, we used mathematical modelling to predict plasticiser loss inside the stagnant environment of a closed container, in the absence of any airflow, at mild temperatures. To the best of our knowledge, plasticiser loss has not been investigated before at mild temperatures inside an enclosure, a scenario that mimics the effect of a closed environment on the degradation of plastic artefacts in storage in museums or more generally on thin plastic samples inside closed containers.
Theoretical plasticiser loss data from the model are compared to Fourier Transform Infrared (FTIR) and 1 H Nuclear Magnetic Resonance (NMR) spectroscopic data from accelerated ageing experiments inside an enclosure at 70 • C and 50% relative humidity (RH) to demonstrate the appropriateness of the model for describing enclosed systems. This is followed by a selection of appropriate ranges in which to vary the diffusion coefficient of plasticisers in different plastics, D s , and their partition coefficients, K, used to simulate plasticiser loss under theoretical scenarios from samples of different thickness, h s , and in enclosures of different size, h g . This analysis investigates the role of the key model parameters on plasticiser loss. Finally, the usefulness of the model in a heritage context is discussed.

Mathematical model
The model describes plasticiser loss from a polymeric slab of uniform (initial) plasticiser concentration placed inside an enclosure. The characteristic dimensions of the slab are L (length), W (width) and h s (thickness), with h s /L≪1 and h s /W≪1; therefore, the object can be  One dimensional, two-phase transient diffusion of plasticiser from plastic which a) rests upon an impermeable surface inside an enclosure (single-sided diffusion) and b) is placed in the middle of the container (symmetrical diffusion from the mid-plane plane located at z = 0). Tables 1 and 2 summarize key parameters and variables of the system. regarded as thin. For this reason, we considered a one-dimensional analysis of the plastic-air system.
In our analysis, we considered low molecular weight, volatile plasticisers such as DEP, DMP and DBP, i.e., those that evaporate rapidly from the surface, rather than deposit on it (which would be the case for the commonly used higher molecular mass PVC plasticisers, such as DEHP). We assumed that initially (that is, at t = 0) the plastic is homogeneous, thus plasticiser is uniformly distributed within the plastic volume. Our system consists of only two layers or phases: the solid plastic and the air. Plasticiser was assumed to migrate from the plastic through a two-phase transient diffusion process. Specifically, plasticiser diffuses from the solid plastic sample across the interface into the air inside the enclosure ( Fig. 2(a)). We also assumed that the temperature is uniform, so that there are no convective effects on gas transport. Finally, no condensation, adsorption or absorption mechanisms are considered on the container walls or on its lid.
To describe the process, we write the mass balance equations for the plasticizer in the two phases (in one dimension) (See Fig. 2(a)) and apply Fick's law of diffusion; this yields: D s and D g represent the diffusion coefficients of plasticiser in the solid phase and in the air, respectively. C s (z, t) and C g (z, t) are the plasticiser concentrations in the solid and in the air expressed in terms of mass per unit volume, so that they are easily comparable to those obtained experimentally. In addition, h s represents the thickness of the plastic, and h g is the distance from the plastic-air interface to the upper enclosure surface (this represents the size of the enclosure). h tot is the height of the container. Table 1 summarizes all key model parameters. For a thin and flat polymeric slab inside an enclosure, we assumed that one surface is resting upon an impermeable solid surface, so here the mass flux of the plasticiser is set to zero. This results in the following boundary condition:

∂C s ∂z
The other surface of the polymeric slab allows plasticiser to diffuse out of the solid into the air. The mass flux of plasticiser that leaves the solid polymer is equal to the mass flux of the plasticiser that enters the air; therefore, we can write: Also, the plasticiser in the air cannot escape the enclosure, which is assumed to be perfectly sealed. Thus, the mass flux of plasticiser at z = h tot is zero, so that:

∂C g ∂z
There is no barrier to transport across the plastic-air interface and local equilibrium prevails at all times. The surface plasticiser concentration in the solid phase and the surface plasticiser concentration in the air are related through a partition coefficient K [46,47] as follows: We also assumed that at t = 0 the plasticiser concentration in the plastic is uniform, that is, C s (z, t = 0) = C 0 , while the plasticiser concentration in the air is zero, C g (z, t = 0) = 0. We solved the system of Eqs. (1)-(6) by using the gPROMS Proc-essBuilder v 1.5.1 -gPROMS -core 6.0.4 2019- 11-13 (55,142) software, by applying the central finite difference method (CFDM). We applied a uniform grid to the solution domain [0, h tot ]. We checked that our solution was grid-independent by comparing simulations where the solution domain had grids with different number of nodes (from 100 to 1000 nodes). We also checked that the solution was independent of the time step used, by comparing simulations with different time steps (from 10 to 10 7 s).
This model describes plasticiser loss under the scenario of singlesided diffusion, which mimics the loss from thin and flat plastic artefacts placed inside an enclosure with one of their surfaces resting upon an impermeable (flat) surface of the container and the other surface exposed to an air gap. We employed the single-sided diffusion analysis to model DEP plasticiser loss from 20DEP/CA samples placed in closed tubes (Section 2.3.2) for a comparison of theoretical loss from the model with plasticiser loss as measured by ATR-FTIR analysis (Section 3.1.2). It was also used to simulate plasticiser loss in the theoretical scenarios of Section 3.3. In Section 3.1.1, where we compare DEP loss as predicted by the model with DEP loss as measured by ATR-FTIR and 1 H NMR analysis, for 20DEP/CA samples placed in closed bottles (Section 2.3.1), the model considers a symmetrical diffusion process from each surface of the samples ( Fig. 2(b)). This is due to the experimental setup allowing for plasticiser being symmetrically lost from both the upper and lower surfaces of the plastic samples. Due to the symmetry, we considered only half of each sample (hence, a thickness of h s /2), with the boundary condition described in Eq. (3) holding at the sample mid-plane. In this case, h tot represents half of the height of the container. By doing so, the problem could be investigated as a single-sided diffusion problem for half of the sample. The analysis is the same whether we consider the volume between the mid-plane and the upper surface or the volume between the mid-plane and the lower surface of the plastic. These considerations are further discussed in Section 3.1.

Application of the numerical solution in defining plasticiser mass loss and theoretical equilibrium times
The profiles of the plasticiser concentration resulting from the numerical solution of Eqs. (1)-(6), C s (z, t), are used to obtain the total plasticiser mass at time t, M tot (t), present in the sample (for the onesided diffusion case) or half of the sample (for the two-sided diffusion case). This is given by: Table 1 Description of the model parameters.

L (m)
Length of the plastic W (m) Width of the plastic hs (m) Thickness of the plastic hg (m) Distance from the plastic-air interface to the upper enclosure surface htot (m) Height of the container (plastic resting on impermeable surface in the container/single-sided diffusion, Fig. 2(a)) Height of half of the container (plastic placed in the middle of container/symmetrical diffusion from both surfaces, Fig. 2 Diffusion coefficient of plasticiser in the plastic Partition coefficient of the plasticiser concentration between the plastic and the air, Cs/Cg where L and W are the length and the width of the plastic object, respectively, and H is the thickness of interest. For a plastic which sits into a container onto an impermeable surface as in Fig. 2(a), H is equal to the sample thickness h s (closed tubes experiment, Section 2.3.2, and theoretical scenarios of plasticizer loss, Section 3.3). In this case, the total plasticiser mass at time t inside the whole plastic is M tot (t). For a plastic which undergoes symmetrical plasticiser loss from both sample surfaces as in Fig. 2(b), H is equal to half of the sample thickness, h s /2 (closed bottles experiment, Section 2.3.1). In this case, M tot (t) represents the plasticiser mass present at time t in half of the sample. We also defined as total mass loss at time t, T L (t) (%), the difference between the total plasticiser mass in the plastic at the time of interest t, M tot (t), and the initial total plasticiser mass, M tot (t 0 ), divided by the initial total plasticiser mass: Table 2 summarizes the defined plasticiser loss terms. The total mass loss at time t, T L (t), is used to describe plasticiser loss in the theoretical scenarios in Section 3.3.
Regarding the theoretical equilibrium times for physical systems controlled by diffusion, which were also used in the analysis of the theoretical scenarios, the transition of the system to the equilibrium state requires mathematically an infinite time [48,49]. But there is a finite time at which the physical system has practically reached equilibrium [50]. These finite times, in a one-dimensional system, are proportional to h 2 /D, where D is the diffusion coefficient and h is the thickness of the medium [50][51][52]. This analysis comes with limitations, because it applies only for homogeneous media with a constant D [51]. In heterogeneous media, where the diffusion coefficient D changes across the different layers, the transition times characterizing the entire physical system are more complex [53]. Thus, we have investigated whether those transition times are valid for our heterogeneous system. However, our system is described by different boundary conditions to those described in [53], so those formulas do not hold here. Therefore, in our analysis we used the simple relation, τ (s) = H 2 /2D s [51].

Table 2
Description of the model variables and other terms.

Symbol (units) Description
Equation

Concentration of plasticiser in the plastic
Concentration of plasticiser in the air Total plasticiser mass in the whole sample (H = hs) or half sample (H = hs/2) at time t

Artificial ageing of cellulose acetate samples plasticised with DEP inside a Duran bottle
Sample preparation. Cellulose acetate (CA) samples were prepared by dissolving 96 g of commercial CA (Aldrich) with an average degree of substitution (DS) equal to 2.48 in 350 mL of acetone during reflux for 4.5 h. Afterwards, with stirring continuing for the entire process, the solution was cooled for 1 h and plasticiser was added. Samples containing 20 wt.% diethyl phthalate (DEP) (99.5% purity, purchased from Sigma Aldrich) were prepared. After 30 min of stirring, the plasticised mixture was poured over a glass dish 26 cm in diameter and covered with a glass lid. After a week of slow evaporation at room temperature, the CA was placed in a vacuum oven to dry for 24 h. Afterwards, the plastic was cut in 2 cm x 2 cm strips with a thickness of 2.2 mm and again placed in the vacuum oven for another 48 h to eliminate acetone. All the samples were stored in the fridge at 5 • C. Samples were named XY/CA, where X refers to the initial amount of plasticiser (in % wt.) and Y to the type of plasticiser (DEP) e.g. 20DEP/CA.
Artificial ageing. Each sample was placed in a 100 mL Duran bottle hanging from the top of the bottle on a woven stainless-steel Grade 304 wire cloth 10 mesh with 1.98 mm aperture (Fig. 3). Each sample aged using a saturated solution of magnesium chloride (MgCl 2 ), sodium bromide (NaBr) or potassium chloride (KCl) (Fig. 3 ATR-FTIR spectroscopy. A Bruker Alpha spectrometer with a Platinum ATR single reflection diamond as Internal Reflection Element (IRE) accessory attached was used. The analysis was performed over a wavenumber range of 400-4000 cm − 1 , using a wavenumber resolution of 4 cm − 1 and 32 scans and the spectra of the samples were collected in absorbance mode. A background spectrum was collected under the same conditions before analysis. Five spectra on the same sample were performed and the resulting spectra were processed using OPUS 7 software.
Plasticiser concentration at the surface. The plasticiser concentration at the sample surface was determined using FTIR spectroscopy through the plasticiser-to-CA intensity ratio (I PL /I CA ). The CH benzene band at 748 cm − 1 was selected to monitor DEP and it was normalised against the peak at 602 cm − 1 , corresponding to the C-C-C backbone of the cellulose ring [14,55,56]. This peak was considered as an internal standard, a constant that is unaffected by the degree of plasticiser and the hydrolysis process [16]. The baseline of the spectrum was defined by averaging the intensity of the spectrum between 2000 and 2200 cm − 1 as this region does not contain any relevant peaks. This value was then subtracted from the intensity at every wavenumber.
I DEP refers to the integrated area of the DEP methyl triplet which is calculated between 1.40 and 1.15 ppm for samples in which DEP was used as the plasticiser, while I IS refers to the integrated area of the Internal Standard (IS) singlet resonance, calculated between 8.57 and 8.37 ppm; N represents the number of 1 H nuclei that correspond to those peaks; M is the molecular mass in g⋅mol − 1 of the compound; m represent the mass of sample (S) and internal standard (IS) used in the analysis and, finally, P IS denotes the IS mass fraction purity which in this work is 99.82%. N for the DEP and IS molecules are equal to 6 and 1 respectively.

Artificial ageing of CA samples plasticised with DEP inside closed aluminium tubes
Sample preparation. Sample preparation was the same as in 2.3.1 for CA samples plasticised with 20 wt. % DEP. The drying procedure took place in a vacuum oven (150 mbar) at 20 • C for 24 h followed by 25 • C for a further 24 h. The samples were then placed in a desiccator until they were cut into 1 × 1 cm 2 pieces of 2 mm thickness 24 h prior to aging. For the rest of our analysis, we will refer to them as closed tubes (at the bottom of the Binder chamber, (Fig. 5)). Other tubes with varying degrees of ventilation were also used but will not be discussed in this paper.
Inside the Binder KBF115 environmental chamber (Fig. 5), a constant climate was achieved by use of an always-on fan (fan speed = 100%) and continuous temperature and relative humidity control. Two iButton DS1923 Hygrochron Temperature/Humidity Loggers were placed, one inside the middle of the closed tubes and one inside the chamber to confirm the desired conditions of 70 • C and 50% RH throughout the experiment. The samples appeared visibly unchanged after aging.
ATR-FTIR spectroscopy -plasticiser concentration at the polymer surface. A Bruker Alpha spectrometer, with a Platinum ATR single reflection diamond as Internal Reflection Element (IRE) accessory attached, was also used for this experiment. The analysis of the received spectra is the same

Applicability of the model to experimental data
We applied the Model Validation tool from gPROMS ProcessBuilder software to fit the numerical values a) of the plasticiser concentration at the plastic surface obtained from the simulations to the experimental values obtained from ATR-FTIR and b) of the total plasticiser concentration in the whole plastic obtained from the simulations to the experimental values found from 1 H NMR data. The fitting yielded the best estimates for the diffusion coefficient D s of DEP in CA samples and for the partition coefficient K at the polymer-air interface of 20 DEP/CA samples. We calculated the value of D g at 70 • C using its value at 25 • C (reported in the literature) and the temperature dependence suggested by Fuller [57,58], D(298) = D(T) 298 1.75 T − 1.75 . Since the measured value at 25 • C is 5 × 10 − 6 m 2 s − 1 [57,59], we obtained a value at 70 • C of 6.4 × 10 − 6 m 2 s − 1 .
The best fit was chosen by considering the Goodness of the fit test of the Model Validation tool of gPROMS ProcessBuilder, a chi-squared test comparing the weighted residual (that is, the difference between the experimental value of the dependent variable and the value predicted by the model) and the expected weighted residual, and by selecting the fitting with the lowest chi-squared value.

Applicability of the model to experimental data
In this section, we applied the model to fit the amount of plasticiser in CA samples plasticised with DEP to experimental values of DEP plasticiser as calculated from ATR-FTIR and 1 H NMR spectra from the artificial ageing experiments described in Section 2.3. The purpose of this analysis is to illustrate that a one-dimensional model of a two-phase transient diffusion process used with appropriate values of the parameters D s , D g and K can successfully describe plasticiser loss from thin plastic samples stored inside closed environments.
For this purpose, we fitted the values of DEP concentration at the surface of the plastic to the experimental surface values as calculated from ATR-FTIR spectra and the numerical values of the total DEP concentration in the whole CA sample to the experimental bulk values of DEP as calculated from the 1 H NMR spectra. The fitting allowed estimating D s and K (for which data are lacking in the literature), while D g was fixed to the value of 6.4 × 10 − 6 m 2 s − 1 .
Given the dimensions of the samples used in the experiments (2 cm x 2 cm x 2.2 mm and 1 cm x 1 cm x 2 mm), we assumed that all plasticiser is lost in one direction, from the largest surface of the sample, as described in Section 2.1 ( Fig. 2(a) and (b)), with loss from the sides of the samples being negligible. Furthermore, while in the gas phase diffusion of DEP occurs in three dimensions in the experiments, given that in the gas phase the diffusion coefficient of DEP is several orders of magnitude larger than in the solid phase, we can expect that on the time scales of interest the gas phase in the bottles or tubes is nearly fully mixed in all three dimensions. Thus, in both phases a one-dimensional analysis is justified.
Moreover, in the closed-bottle experiments as well as in the tube experiments, owing to the system geometry, the plasticiser could leave the sample from both the upper and lower surfaces. In the former case, the sample was (essentially) placed in the middle of the bottle, so there were two equal volumes of gas above and below the object. Consequently, the diffusion problem was symmetrical with respect to the middle plane of the sample (i.e., the plane located at the same distance  from the top and bottom surfaces of the bottle, Fig. 2(b)). Owing to symmetry, the diffusive flux at the mid-plane of the sample must be zero, and therefore the boundary condition in Eq. (3) is valid. Furthermore, Eq. (4) holds both at the upper and lower surfaces. Under these assumptions, we simulated only half of the sample (the part between the mid-plane and one of the two surfaces), considering only half of the sample thickness (i.e., h s /2). Accordingly, as discussed in Section 2.2, in these simulations M tot (t) represents the mass of plasticiser contained only in half of the sample. Conversely, in the closed-tube experiments, some plasticiser is lost from the bottom surface of the sample, for this is not adjacent to the tube wall, but the volume of gas between the bottom surface of the sample and the tube wall is negligible compared with the volume of gas between the top surface of the sample and the tube wall. Hence, neglecting the contribution of the first volume, and employing Eq. (3) as boundary condition at the bottom surface, is acceptable.

The closed-bottle experiment
We first consider the closed-bottle experiment described in 2.3.1. Fig. 6 illustrates the fitting of the numerical values of DEP concentration at the surface of the plastic (black squares) to the experimental values of DEP (red circles) as calculated from ATR-FTIR spectra (as described in Section 2.3.1). The model considers a single-sided diffusion for DEP ( Fig. 2(a)). The numerical values of DEP concentration at the plastic surface were defined as: dM zIR pen (t) dM zIR pen (t) + dm CA, zIR pen * 100 (wt.%) where z IR pen is the z value inside the 20DEP/CA samples at a distance d IR from the plastic-air interface that corresponds to the penetration depth of IR radiation in CA. This is about 3 µm, measured from the plastic-air interface, thus d IR ∼3 µm [60]. Additionally, dM zIR pen (t) is the DEP mass in the CA sample within the IR radiation penetration depth; thus, inside the volume dV = LWd IR :

dM zIR pen (t) = C s, zIR pen (t)dV
where C s, zIR pen is the plasticiser concentration in the plastic at z IR pen . We assumed that it is uniform inside the volume dV. Furthermore, dm CA, zIR pen is the CA mass within the IR radiation penetration depth, thus inside the volume dV. The initial DEP concentration, in wt.%, is 20%. We denoted this as g 0 . Hence, g 0 = 0.2. It was also assumed that the mass of CA is constant; thus, the total CA mass inside the sample is m CA = (1 − g 0 ) m 0 , where m 0 is the initial total sample mass. Then, the CA mass inside the volume dV is:   (11) where V sample = LWh s is the sample volume. The fitting yielded the value D s =4.1 × 10 − 11 m 2 s − 1 and K=(1.440±0.008)x10 2 . However, the equilibration of the plasticiser concentration at the surface is very rapid (after 1 d) and the estimated D s value through the ATR-FTIR data in Fig. 6 derives from two data points only (day 0 and day 1); accordingly, we regarded this value as unreliable and for the rest of our analysis we did not employ it. In Fig. 7 the numerical values of the total DEP concentration in the whole plastic from simulations (black squares) were fitted to the experimental values of total DEP (red circles) as calculated from 1 H NMR data by using Eq. (9). The model here considers a symmetrical diffusion of DEP from the 20DEP/CA samples (Fig. 2(b)). D g was fixed to the literature value of 6.4 × 10 − 6 m 2 s − 1 , as described earlier. K was also fixed to the value of 10 2 , similar to the value deriving from the previous fitting by using the surface ATR-FTIR data. The numerical values of the total DEP concentration in the whole plastic, DEP tot were defined as: Due to the symmetry of the system, 2M tot (t) is the total plasticiser mass in the whole polymer.
The fitting yielded the value D s =(1.510±0.027)x10 − 13 m 2 s − 1 . We also investigated a fitting of the 1 H NMR data by considering as fitting parameters both D s and K. The resulting values are very similar: D s =(1.312±0.035)x10 − 13 m 2 s − 1 and K=90.0±1.80.

The closed-tube experiment
We also fitted the numerical values of DEP concentration at the surface of the plastic sample obtained from the simulations to the experimental values calculated by using the plasticiser-to-CA intensity ratio (I PL /I CA ) from the ATR-FTIR analysis from the closed-tube experiment (Fig. 8). The numerical values at the polymer surface were calculated according to Eq. (10). D g was fixed to the literature value of 6.4 × 10 − 6 m 2 s − 1 . D s was fixed to the value of 1.51 × 10 − 13 m 2 s − 1 deriving from the fitting from the closed-bottle experiment, since the temperature and relative humidity were the same and the samples used in each experiment were produced in the same way. The fitting provided K = (2.486±0.061)x10 3 . Here, we also explored a fitting of the ATR-FTIR data by considering as fitting parameters both D s and K. The fitting yielded the similar values D s = (2.471±0.032)x10 − 13 m 2 s − 1 and K = (2.013±0.045)x10 3 .

Discussion on the values estimated from the fitting
In the closed-bottle experiment, we estimated D s from fitting the numerical values of total DEP to 1 H NMR data. The fitting provided D s = 1.51 × 10 − 13 m 2 s − 1 . The estimation of D s from 1 H NMR data is regarded as reliable, as it is based on several experimental data points and the 1 H NMR analysis is applied to the whole sample.
Furthermore, the ATR-FTIR data are useful for estimating K, which relates the plasticiser concentrations in the two phases at the sample interface. The estimated value of K ≈ 10 2 from fitting the numerical DEP concentration values at the polymer surface to the experimental ones from ATR-FTIR data also works well when we fitted the DEP concentration in the whole plastic obtained numerically to the experimental DEP concentration from 1 H NMR data.
However, the estimated value K = (2.486±0.061)x10 3 obtained from the closed-tube experiment is higher compared to the value K ≈ 10 2 obtained from the closed-bottle experiment. This could be attributed to the different method used to produce the environment inside the bottles. The salt solution used in the bottles to control the relative humidity could act as a sink, absorbing plasticiser from the air. Indeed, recent experimental work in our group has identified the presence of DEP in salt solutions employed for artificial ageing experiments using UV-Vis spectroscopy, demonstrating that the salt solutions act as a sink for DEP (unpublished research from our group). This phenomenon would result in a higher plasticiser loss, resulting in a lower K value compared to that obtained from the closed-tube experiment.
Overall, by fixing D g at 70 • C according to the literature, we have managed to successfully fit the numerical values of DEP concentration at the polymer surface to the experimental values calculated from ATR-FTIR data and numerical values of total DEP concentration in the whole plastic to experimental values of DEP obtained from 1 H NMR data by using similar K values (K = 1.44 × 10 2 and K = 10 2 ). The K values estimated in the current work should be considered as estimates of the partition coefficient of DEP between the plastic samples and the surrounding environment. We have also fitted the DEP concentration values from the simulations to the experimental values from distinct experiments conducted at the same conditions of temperature and relative humidity by using the same value for the diffusion coefficient D s .
Thus, a one-dimensional model for two-phase transient diffusion, if used with appropriate values for the parameters D s , D g and K, can successfully describe plasticiser loss from thin plastic samples stored in closed environments. We discuss the appropriateness of the estimated values in Section 3.2, where we compare the estimated values with the real values of D s and K from related systems.

Selection of parameter ranges
In this section, we identify appropriate ranges for the model parameters D s , D g , K, h s and h g that will be used for future analysis.
The diffusion coefficients of additive-like molecules of various molecular weights in different thermoplastic materials in both rubbery and glassy states are seen in Fig. 9 [61]. We observe that for low molecular weight plasticisers like DEP, dimethyl phthalate (DMP) or dibutyl phthalate (DBP) with a molecular weight ~ (200-300) g/mol, their diffusion coefficients are within a range 10 − 21 m 2 s − 1 < D s < 10 − 14 m 2 s − 1 in the glassy state. We used the lower and upper boundaries of this range to predict total plasticiser loss in the theoretical scenarios in Section 3.3.
The estimated value of D s from the experiments, D s = (1.510±0.027) x10 − 13 m 2 s − 1 , is not within these boundaries. Having checked that 20DEP/CA samples after the artificial ageing of 2.3.1 were still in the glassy state, we attribute the higher value estimated for D s to the higher temperature of 70 • C at which we conducted our experiments and to the different plastic system compared to those to which Fig. 9 refers. Diffusion coefficients are known to increase with temperature [62].
In the following theoretical scenarios, we set D g = 5 × 10 − 6 m 2 s − 1 . This value refers to the diffusion of DEP in air in mild environments at 25 • C, as an example of loss of a low molecular weight plasticiser from a polymer in realistic environmental conditions. The range of 10 − 6 <D g <10 − 5 m 2 s − 1 was identified as representing the values of the diffusion coefficients in air for various plasticisers more broadly. Note that since D g ≫D s , the value of D g used in the model, if taken within this range, does not affect the behaviour of our system significantly. This is also illustrated in detail in Section 3.3.
Regarding the partition coefficient K, this dictates the ratio of the equilibrium concentrations of plasticiser between the plastic and gas phases at the sample interface and is dimensionless (if the concentrations in the two phases are expressed in the same units). We selected a range of K values between 10 3 and 10 7 , which have been measured for volatile organic compounds of various molecular weight from vinyl flooring at 20 • C [63]. We used this range to investigate the effect of K on total plasticiser loss. Furthermore, the lower and upper boundaries of this range were used to predict the total plasticiser loss in the theoretical scenarios considered in Section 3.3.
The estimated K value of 10 2 obtained from the closed-bottle experiment described in 2.3.1 lies outside this range. A Van't Hoff dependence is known for K [64], so that K increases with the temperature. However, the system inside the bottles is quite complex, since it includes salt solutions to control the RH. This system may act as a plasticiser sink, absorbing plasticiser from the polymer surface, a phenomenon that the model neglects. As also explained above, experimental work in our group has identified the presence of DEP in salt solutions used for artificial ageing experiments using UV-Vis spectroscopy. This shows that the salt solutions act as a sink for DEP. In addition, plasticiser may be adsorbed in the plastic cap or on the glass wall [65], which again is not considered in this model. Conversely, the estimated K value of (2.486±0.061)x10 3 from the closed-tube experiment described in 2.3.2 is within the literature range. This system is simpler and does not include the salt solution inside the tubes to control the relative humidity while the cap and the wall of the tube are of the same material.
We also assumed that the thickness of the thin polymeric plasticised sample varies between 10 − 3 <h s <10 − 2 m and the distance of solid-gas interface from the upper surface of the enclosure varies between 10 − 2 <h g <10 − 1 m. In a museum conservation context, discussions with heritage professionals and surveys of plastic artefacts (unpublished work undertaken within our research group) have shown that plastic artefacts plasticised with DEP that exist in museum collections are mainly thin, of a few mm thickness. Moreover, they are mainly preserved inside enclosures where the distance of the upper part of the artefact from the lid of the container is either a few mm or a few cm, without exceeding 10 cm.

Theoretical scenarios of plasticiser loss
Using the values discussed above, we simulated the plasticiser concentration over time across the samples with different thicknesses inside differently sized enclosures for single-sided diffusion ( Fig. 2(a)). The results are shown in Table S1. Fig. 10 reports the reduction of the plasticiser concentration (i) at the mid-plane (0.5 mm) of the sample, (ii) at a distance of 2 µm from the solid-gas interface, and (iii) at the interface of a 1 mm sample stored inside a 10 cm enclosure in the scenario where D s = 10 − 14 m 2 s − 1 and K = 10 3 for short ( Fig. 10(a)) and long ( Fig. 10(b)) times. The initial plasticiser concentration is taken to be 20 wt% as in the closed-bottle experiment.
The dip of the black dot curve representing the surface plasticiser concentration at short times ( Fig. 10(a)) is attributed to the initial conditions of the system, where the initial concentration in the sample (C 0 ) is 292.63 kg m − 3 , while the initial plasticiser concentration in the air is zero. Equilibration at early times results in a rapid initial surface loss, the interface concentration in the solid sample reaching ∼44 kg m − 3 after 1 min. The initial rapid depletion at the surface is caused by surface plasticiser loss being faster than diffusion in the bulk of the solid. This is also illustrated in Fig. 11, which shows the plasticiser concentration profile across the 1 mm sample after 1 d, 5 d, 15 d, 30 d and 1 yr inside the 10 cm enclosure, for the same D s = 10 − 14 m 2 s − 1 and K = 10 3 . The high concentration gradient for z ∈ (0.9, 1) mm at 1 d (blue cross) is related to the high loss from the surface at short times, in contrast to the bulk concentration which has not yet been affected.
However, the concentration at the plastic surface does not remain at such low values in a closed system. As time evolves, surface loss is replenished by diffusion from the bulk of the solid. So, the plasticiser concentration in the bulk reduces. As plasticiser accumulates in the gas phase inside the enclosure and the system approaches equilibrium, the loss rate from the surface reduces. Therefore, the concentration at the surface increases over time, while inside the bulk it decreases (Figs. 12, 5d, 15d, 30d). Thus, the plasticiser gradient at the surface decreases over time.
As the system reaches equilibrium in the gas phase, the concentration in the polymer levels out at a concentration lower than the initial value but larger than the minimum value attained for short times (Fig. 10(a)). Table S1 shows the outputs of a range of simulations using the parameter values discussed in Section 3.2 and the outputs described in Table 2, together with the theoretical equilibrium times τ (see above).
We also assume that all the plastic samples are new; thus, they have an initial uniform plasticiser concentration C 0 = 292.63 kg m − 3 , a value that corresponds to 20 wt% DEP in the CA samples used in our experiments.
Despite referring to DEP, this analysis could also represent other plasticisers with MW ~ (200-300 g/mol), since the range of D s values used is broad 10 − 21 <D s <10 − 14 m 2 s − 1 (see Section 3.2).
We have also observed non-uniform plasticiser concentration across cross-sections of aged CA samples, as measured by using Fourier Transform Infrared (FTIR) microscopy (Fig. 12) in a recent publication from our research group [66]. The CA samples in Fig. 12 were aged for 1 or 3 days at 70 • C at different RH conditions (50% and 80%). The initial plasticiser concentration across the samples was 20 wt.%. For all the RH values, we observe a plasticiser concentration profile across all the aged samples, with lower concentrations in the layer between 0 and 10 µm (8-17 wt.%) and higher ones in the layer between 10 and 200 µm (13-20 wt.%). In the layer between 50 and 200 µm, the plasticiser concentration remains almost constant at higher values (15, 18 and 20 wt.%). This observation suggests that plasticiser loss from 20DEP/CA samples is diffusion-controlled at temperatures of 70 • C. Thus, the profiles we simulated in Fig. 11 highlighting a diffusion-controlled system are in good agreement with the ones we measured experimentally by using Fourier Transform Infrared (FTIR) microscopy. In the same work, we observed that plasticiser loss is also diffusion-controlled for historical CA samples which have naturally aged (Fig. 4 in reference [66]).
The analysis referred to in Table S1 explores the role and importance of the model parameters K, D s , h s and h g on plasticiser loss, under the specified ranges of values. By comparing the total loss between the different scenarios, we conclude that the partition coefficient K plays a key role in controlling total loss. When K rises from 10 3 to 10 7 , so by a factor of 10 4 , the total loss decreases by approximately the same ratio ∼10 4 (compare Sim 1-2, 5-6, [9][10][13][14]. This is also shown in Fig. 13, where the plasticiser concentration at the mid-plane of the 1 mm thick polymer with D s = 10 − 14 m 2 s − 1 is reduced by 1% for K = 10 4 but by 10% for K = 10 3 . The equilibrium times are similar for the investigated K values. Diffusion coefficients do not affect the total loss in this system, when the other parameters do not change. When D s decreases from 10 − 14 to 10 − 21 m 2 s − 1 , the total loss remains the same (compare Sim 1-3, 2-4, 5-7, 6-8, 9-11, 10-12, 13-15, 14-16). D s controls the equilibrium times, or how fast the loss process is. This is illustrated in Fig. 14 for the 1 mm plastic sample stored inside the 10 cm enclosure, for K = 10 3 . While D s increases from 10 − 15 m 2 s − 1 to 10 − 14 m 2 s − 1 , the equilibrium time decreases from 10 yr to 1 yr, while the concentration reduction is the same. This analysis agrees with the estimation of equilibrium times according to the approximation τ = h 2 s /2D s , in which the equilibrium time depends only on D s if the sample thickness is the same.
The thickness of the plastic and the enclosure size affect the total loss. This is proportional to the enclosure size under the investigated range of values and between the same K, D s and h s scenarios (compare Sim 1-5, 2-6, 3-7, 4-8). It is also inversely proportional to the plastic thickness under this range, for the scenarios of the same h g , K and D s (compare Sim 1-9, 2-10, 3-11, 4-12, 5-13, 5-14, 7-15, 8-16).
Additionally, the significance of the thickness of the plastic is highlighted by the equilibrium times, which agree with the expression τ=h 2 s /2D s . As it is illustrated in Fig. 15, the equilibrium time for a 1 mm plastic is ∼1.7 yr, while for a 1 cm plastic is ∼160 yr inside the 10 cm enclosure when K = 10 3 and D s = 10 − 14 m 2 s − 1 . The expression τ=h 2 s /2D s yields 1.6 yr for the 1mm plastic and 160 yr for the 1cm plastic, values that are close to those estimated from the simulations.
Finally, as Fig. 16 shows, D g does not affect the plasticiser concentration over time within the examined range of values. Thus, in the range of values considered for D g , the total loss does not depend on its value.

Discussion on the usefulness of the model
Our model offers insight into the degradation of thin plastics due to plasticiser loss inside an enclosure. The work focuses on the migration of low molecular weight plasticisers from thin polymeric materials of different compositions with different storage dimensions and material thicknesses. The model describes how the system evolves in time.
In a conservation context, storing thin plastic artefacts inside small containers is a common method used by museum professionals to limit degradation. Our model allows them to estimate the lifetime of plastic artefacts under different storage scenarios. This analysis can inform conservation guidelines about storage dimensions and strategies. It could also prompt discussions inside the organisation regarding the effectiveness of current storage dimensions and generate options for improving storage. This could spark further research and development of innovative methods to minimize plasticiser loss.
Similarly, the time evolution of plasticiser loss, accompanied by the knowledge of equilibrium times, could facilitate decisions about how to conduct condition assessments, or shipping, or exhibiting the artefact. If equilibrium has already been established, then the further plasticiser loss resulting from opening the container, and the consequent risk of damage and loss of value, especially for the surface layers of the object, should be taken into account. This issue is relevant when considering moving objects from one container to another. Additionally, the loss of information, from not opening the container, and the loss of historical and cultural representation, from not displaying the artefact, should be considered. This highlights the contribution that our model may have to risk assessment analysis. Table S1, which summarizes total loss and equilibrium times, could be useful if correlated with the risk of damage for each scenario and loss of value. On that account, the model could assist decision making in conservation and inform conservation planning.
Moreover, artefacts that have been treated on their surfaces may have undergone compositional alterations there. This may change the partition coefficient value K, allowing for more or less plasticiser to escape from the surface, depending on the treatment method. Our analysis in this paper has also underlined the importance of the K value in terms of plasticiser loss. This work has also identified gaps in the literature regarding K values for plasticisers and plastic/air interfaces, highlighting the need for further research on this subject. Finally, our approach has possible implications to packaging design regarding thin plastic products with low molecular weight plasticisers, due to enabling predicting plasticiser loss for different enclosure sizes and plastic systems.

Conclusions
We have developed a model that predicts plasticiser loss over time from thin plastic objects stored inside closed environments through a two-phase transient diffusion mechanism. Our model provides insight into the relationship between plasticiser loss rate, equilibrium time, enclosure dimensions and material thickness in various polymers. The model predictions were compared to FTIR and 1 H NMR spectroscopic data from degradation experiments at 70 • C and 50% RH, showing that this modelling approach effectively describes the process, provided that appropriate parameter values are used.
Our analysis has explored the role of K, D s and h g on total plasticiser loss through theoretical scenarios. K controls total loss and D s controls the equilibrium times; however, being material properties, they cannot be adjusted after manufacture, although they change with temperature. So, regarding museum applications, temperature control is crucial, because it affect to the total plasticiser loss at equilibrium and the time that the system takes to reach equilibrium.
Conversely, the dimensions of the storage system, which are also important for plasticiser loss, can be adjusted to inhibit loss. This could spark innovative planning in terms of the dimensions of enclosures for plastic objects that could be designed to minimise plasticiser loss. Ultimately, this model can be applied in museum conservation contexts, or other contexts where thin plastics with low molecular weight plasticisers are stored in closed environments, to predict lifetimes of plastic samples. It can therefore be a useful tool to explore the longevity of thin plastics under different scenarios and facilitate decisions regarding packaging strategies for products or preservation methods and storage systems for plastic museum artefacts. In that respect, the application of the model could contribute to potential inspiration for future conservation guidelines, strategies and policies.

Funding
This project has received funding from the European Research