Chemo-mechanical model for degradation of oil paintings by amorphous and crystalline metal soaps

Metal soap formation is recognised as a critical degradation mechanism in historical oil paintings, which threatens the preservation of museum collections worldwide. Metal soaps form via a complex sequence of chemical reactions between metal ions released by the pigments and saturated fatty acids originating from the drying oil. The latest advances in chemistry research suggest that metal ions and saturated fatty acids may initially react by means of a reversible reaction, which leads to the formation of metal soaps in an amorphous state. Metal soaps may subsequently crystallise via an irreversible reaction into large aggregates that deform the paint layers, potentially triggering delamination, cracking, and ultimately flaking of the paint. This paper proposes a chemo-mechanical model to predict metal soap formation and the consequent mechanical damage in historical oil paintings. The chemical process is described in terms of a set of diffusion–reaction equations, which account for both the reversible reaction between free saturated fatty acids and metal ions forming amorphous metal soap, and the subsequent irreversible reaction to crystalline metal soap. The chemical model is two-way coupled with a mechanical model that effectively describes the cracking processes caused by metal soap formation and growth. The coupling is generated from the mechanical model by accounting for the development of a chemically-induced growth strain in the crystalline metal soap. In addition, the presence of cracks locally hampers the diffusion of chemical species, which is taken into account in the chemical model through a dependency of the diffusion parameter at the crack faces on the amount of mechanical damage generated. The spatial development of the crystalline metal soap phase is simulated by using a tailor-made scanning algorithm that identifies the reaction zone in which metal soap formation takes place. The proposed model is calibrated on experimental data presented in the literature. The model is subsequently applied to analyse two numerical examples that are representative of typical metal soap-related degradation processes observed in historical oil paintings, revealing that the growth process of crystalline metal soap, the deformation of the paint surface, and the consequent cracking and delamination patterns are predicted in a realistic fashion.


Introduction
The process of metal soap formation and growth has been identified by cultural heritage conservators as one of the most crucial degradation mechanisms for the appearance, integrity and longevity of historic oil paintings (Cotte et al., 2017;Noble, 2019).Metal soaps typically appear as opalescent aggregates that deform the pictorial layers, possibly protruding through the paint surface (Higgitt et al., 2003;Noble and Boon, 2007;Cotte et al., 2007;Keune and Boon, 2007).They may additionally lead to an increased transparency of the paint, thereby revealing the colour of the underlying paint layers (Noble and Boon, 2007).Metal soaps do not only influence the visual appearance of oil paintings, but may also promote cracking, delamination and ultimately In the bottom orange-red paint layer, which contains lead white and red lead, three lead soap protrusions have formed.An interface undulation can be observed in this cross-section, whereby the crystals in the bottom layer start to penetrate the layers above (Thoury et al., 2019).(b) Johannes Vermeer, View of Delft, ca.1660-1661, Oil on canvas, Mauritshuis, The Hague, The Netherlands.The detail shows a bright-field image of a paint cross-section taken at the left edge of the painting, where the green (now blue) bush is painted over the red-tiled roof, showing the formation of large metal soap protrusions in the dark-red paint layer that have erupted through the blue-green surface paint layer ( van Loon, 2009; van Loon et al., 2012).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Oil paintings consist of layers of a drying oil binding medium (typically linseed oil), which contain pigment particles and adhere to a substrate made of prepared canvas or wood.From a chemical point of view, drying oils are triglycerides, composed of three fatty acid molecules bound to a glycerol backbone via an ester bond.Fatty acids are carboxylic acids, i.e., organic compounds of carbon, oxygen and hydrogen, with long hydrocarbon chains (Lazzari and Chiantore, 1999).Most of the fatty acids chains in a drying oil are unsaturated, containing one or more C=C bonds.However, saturated fatty acids, such as palmitic and stearic acids, are also present in the oil binder, representing approximately 5%-15% of the total fatty acid content.During paint drying, the unsaturated fatty acid chains in the oil binder react by means of autoxidation reactions to form a densely cross-linked polymer network (van Gorkum and Bouwman, 2005;Baij et al., 2019;Bonaduce et al., 2019).Moreover, on a larger time scale, due to the presence of water molecules in the paint film, the ester bonds between the fatty acid chains and the glycerol backbone can be broken by means of hydrolysis reactions (Baij et al., 2019;van den Berg et al., 2001;van den Berg, 2002;Bonaduce et al., 2012).This leads to the release of free saturated fatty acids as a reaction product (Erhard et al., 2005).The process of metal soap formation is especially relevant in the case of lead-or zinc-based pigments, which are present, for instance, in the common lead white and zinc white paints (Keune et al., 2008).Lead and zinc ions that are released from the pigments may become complexed with the carboxyl groups (COOH) of the polymer network, forming an ionomer (Hermans et al., 2015(Hermans et al., , 2016a;;Baij et al., 2018a;Beerse et al., 2020), see Fig. 2(a).The free saturated fatty acids present in the system may further react with the metal ions from the ionomer, forming complexes of amorphous metal soaps, see Fig. 2(b).The formation of amorphous metal soaps is expected to be a reversible process, for which the metal ions can be bound either to the carboxylic acid groups of the ionomer or to those of the free saturated fatty acids (Hermans et al., 2015).If the local concentration of amorphous metal soap becomes sufficiently high, it can subsequently crystallise by means of an irreversible chemical reaction and coalesce into larger metal soap aggregates, reaching diameters varying between 50 and 500 μm (Noble and Boon, 2007), see Fig. 2(c).
Fig. 3 shows the internal material structure of a metal (zinc) soap nucleus at the nano-scale.The image, which has been taken from Hermans et al. (2018), has been obtained by using transmission electron microscopy, and clearly illustrates that the crystalline metal soap dendrites (indicated in red) that have developed inside the nucleus are surrounded by paint material (indicated in blue) based on a linseed oil polymer.In other words, the metal soap nuclei, and also the larger aggregates developing from these nuclei, are not fully occupied by crystalline metal soap, but also contain paint material, as confirmed from other experimental studies (Henderson et al., 2019;Ortiz Miranda et al., 2020).For simplicity reasons, however, in the present work the metal soap nuclei and aggregates will be designated as ''metal soap crystals'', and in the model formulation are effectively treated as homogeneous and isotropic.The latter assumption is reasonable, as the model refers to the length scale associated to the thickness of a paint layer, i.e., the micro-scale, which is considerably larger than the nano-scale considered in Fig. 3. Nevertheless, in the numerical simulations the initial concentration of metal ions will be determined by considering that the aggregates are not fully occupied by crystalline metal soap, which essentially means that only a certain fraction of the metal ions present in the initial paint material is available for metal soap formation.
The process of metal soap formation has been primarily investigated in the literature from a chemical perspective, by focusing on the insitu identification of metal soap in historic paintings (Noble et al., 2002;van der Weerd et al., 2002;Plater et al., 2003;Keune et al., 2011;Chen-Wiegart et al., 2017), or by exploring the kinetics of the chemical reaction (Cotte et al., 2006;MacDonald et al., 2016;Izzo et al., 2021) and the composition of metal soap aggregates (Hermans et al., 2016b).On the contrary, only a few studies are available that investigate the consequences of metal soap formation on the mechanical response of paint films.In van den Berg (2002), it has been reported that metal soap formation can contribute to an increase of oil paint stiffness and brittleness.Moreover, the presence of metal soap aggregates has been associated to micro-cracks and delamination phenomena (Maines et al., 2011;Helwig et al., 2014;Raven et al., 2019;Osmond et al., 2012), which are triggered by the mechanical strain (and thus stress) generated in the painting system under metal soap growth.Recently, a computational chemo-mechanical model has been proposed that predicts damage development in paint layers as a result of metal soap formation and growth (Eumelen et al., 2019).The model departs from a simplified chemical description in which the intermediate state of amorphous metal soap is neglected and only the irreversible crystallisation process is modelled.Accordingly, the diffusion-reaction model is formulated solely in terms of the concentration of saturated fatty acid and the volume fraction of crystalline metal soap formed.The chemo-mechanical model is implemented within a finite element framework, whereby the chemical and mechanical fields are solved in an incremental-iterative fashion, using a staggered solution scheme.
Departing from the formulation presented in our previous work (Eumelen et al., 2019), in the present contribution the chemo-mechanical model for metal soap formation is substantially enhanced by accounting for both the reversible reaction between free saturated fatty acids and metal ions forming amorphous metal soap, and the irreversible reaction governing the process of metal soap crystallisation.Along similar lines as the one-dimensional formulation proposed in Hermans (2017), these chemical processes are described via a set of (diffusion-)reaction equations in terms of the concentrations of free saturated fatty acids, metal ions, amorphous metal soap and crystalline metal soap.Although, in principle, both the free saturated fatty acids and the metal ions can diffuse in the ionomer (Tudry et al., 2012), metal ion diffusion is ignored here, as this process typically is very slow compared to fatty acid diffusion.As a consequence, the free saturated fatty acids are considered as the sole diffusing species, for which a diffusion-reaction equation is formulated, whereas the time evolution of the concentrations of the metal ions and the amorphous and crystalline metal soaps are governed by reaction equations.Along the lines of Eumelen et al. (2019), the coupling between the diffusion-reaction model and the mechanical analysis is obtained by defining a chemically-induced growth strain, which quantifies the effect of the expansion of the growing metal soap crystal on the stress field generated in the painting.Additionally, the local change in the mechanical properties associated to the formation of crystalline metal soap is computed by a basic phase transformation model, whereby the effective stiffness is obtained as the volume average of the properties of the material phases present in the specific material point.Mechanical equilibrium equations finally complete the definition of the mechanical problem.The paint geometry is discretised into plane-strain continuum elements that represent the two-dimensional bulk responses of the metal soap and paint materials.Further, in order to describe the onset and propagation of discrete cracks at arbitrary locations in the paint domain, cohesive interface elements are placed between all continuum elements modelling the paint configuration, in correspondence with the approach originally proposed in Xu and Needleman (1994).The constitutive behaviour of a crack is defined in accordance with the interface damage model proposed in Cid Alfaro et al. (2009).The presence of cracks locally hampers the diffusion of the free saturated fatty acids, which is accounted for by specifying the constitutive relation between the flux and the concentration of free saturated fatty acids across a crack as a function of the mechanical damage.Consequently, the chemo-mechanical model is two-way coupled.
The enhanced chemo-mechanical model proposed in this work is applied to analyse two different boundary value problems, which are representative of typical metal soap-related degradation mechanisms observed in historical paintings.The first boundary value problem focuses on the interplay between the growth process and coalescence of multiple metal soap crystals and the crack patterns generated in a single-layer paint system.The second boundary value problem refers to a multi-layer paint system composed of an upper pictorial layer and a supporting ground layer, and focuses on the competition between surface cracking in the upper layer and delamination at the interface between the two layers, as induced by the formation and growth of a single metal soap crystal in the ground layer.This paper is organised as follows.Section 2 describes the chemomechanical modelling framework, including the diffusion-reaction equations, the mechanical model and the coupling terms between the two physical processes.In Section 3, the chemical properties that serve as input for the diffusion-reaction model are calibrated based on the experimental results presented in Hermans (2017).Section 4 treats the definition of the geometry, material properties, initial and boundary conditions and the numerical results of the single-layer paint model.In Section 5, the geometrical details, input parameters, initial and boundary conditions and simulation results of the multi-layer paint model are discussed.Finally, the main conclusions of the work are summarised in Section 6.

Model description
The chemo-mechanical model describing the degradation of oil paintings by metal soap formation consists of a chemical model and a mechanical model, which are successively described in Sections 2.1 and 2.2.The numerical solution procedure that is used to solve the coupled chemical and mechanical models is reviewed in Section 2.3.

Chemical model
Departing from the formulation presented in Eumelen et al. (2019), in the present contribution the chemo-mechanical model for metal soap formation is substantially enhanced by accounting for both the reversible reaction forming amorphous metal soap, and the irreversible reaction leading to metal soap crystallisation.The proposed model is based on the following assumptions.First, the initial state of the domain is assumed to be represented by an ionomer characterised by given concentrations of free saturated fatty acids and metal ions.Accordingly, the autoxidation reactions leading to the formation of the ionomer and the reactions associated to the release of free saturated fatty acids by the oil binder and the metal ions from the pigments are not explicitly modelled.Moreover, despite that the paint is a heterogeneous material, composed of the ionomer with embedded pigment particles, it is here treated as an equivalent, homogeneous continuum, characterised by effective chemical and mechanical properties.The paint and the fully crystallised metal soap are the two physical material phases that together form the paint system.Conversely, the free saturated fatty acids, the metal ions, and the amorphous metal soaps only enter the model through their concentrations.It is assumed that metal soap crystallisation develops from pre-existing crystal nuclei, which are small regions in which fully crystallised metal soap is present, see Fig. 3.The chemical reactions leading to metal soap formation are enforced to occur in a (moving) transition region named the ''reaction zone'', which is located at the interface between the metal soap crystal and the adjacent paint material.
The chemical process leading to the formation of crystalline metal soap can be expressed by means of two successive chemical reactions (Hermans, 2017): in which M + Fa denotes a metal ion M reacting with  free saturated fatty acid molecules Fa.The term MFa  indicates the metal soap complex, for which prefixes ''a'' and ''c'' refer to the amorphous and crystalline states of the metal soap compound, respectively.Further,  + and  − are the rate constants for the formation and dissolution of amorphous metal soap, respectively, while  crys is the rate constant of the irreversible crystallisation reaction.In the present work leadbased pigments will be considered, for which  = 2.The chemical reactions described by Eq. ( 1) can be converted to a system of diffusionreaction equations, expressed in terms of the individual species that are involved in the chemical process, i.e., saturated fatty acids, metal ions, amorphous and crystalline metal soaps.However, the diffusion process of saturated fatty acids is characterised by a significantly smaller characteristic time scale than that of the metal ions (Baij et al., 2018a) and of the amorphous or crystalline metal soap compounds.Hence, metal ion diffusion and amorphous metal soap diffusion are ignored, and the free saturated fatty acids are selected as the sole diffusing species.
Based on the above assumptions, the chemical process is defined by the following system of coupled (diffusion-)reaction equations (Hermans, 2017): (2) Here [], with [] ≥ 0, represents the available concentration of the chemical species , with  = Fa, M, aMFa 2 and cMFa 2 .Further,  represents time, and ⋅ denotes the divergence operator, with  being the gradient operator.Finally,  indicates the diffusion coefficient.
Refer first to Eqs.
(2) 1 -(2) 3 that describe the development of the concentration of the species  = M, aMFa 2 and cMFa 2 that react but do not diffuse.The minus sign in front of some of the reaction terms   indicates that, if the value of   is positive, the corresponding species  is consumed during the specific chemical reaction.Define now  +  as the reaction rate of the reactants, which characterises the ''forward reaction'', and  −  as the reaction rate of the reaction products, characterising the ''backward reaction''.Following Oxtoby et al. (2008), the reaction terms  aMFa 2 and  cMFa 2 can be expressed as (3) Note that Eq. ( 1) 2 describes the irreversible chemical reaction from amorphous towards crystalline metal soap.Therefore, in Eq. (3) 2 , the reaction term  cMFa 2 is prescribed to be non-negative,  cMFa 2 ≥ 0.
The specific expressions for the forward and backward reaction rates in Eqs.
(3) can be formulated in terms of the concentrations of the involved species, raised to a power equal to their stoichiometric coefficient (Oxtoby et al., 2008;Zumdahl and Zumdahl, 2017).Considering the chemical reactions, Eq. (1), and recalling again that  = 2, leads to (4) Inserting Eqs.(4) into Eqs.(3) finally specifies the reaction terms  aMFa 2 and  cMFa 2 as (5) which need to be substituted in the system of (diffusion-)reaction equations, Eq. ( 2).Consider now Eq.(2) 4 that describes the diffusion and reaction processes of the free saturated fatty acids [Fa].The reaction term  aMFa 2 present in the right-hand side of the equation is defined by Expression (5) 1 , and introduces a coupling with the reaction equations, Eqs.(2) 1 -(2) 3 , of the non-diffusing species.Note further that the factor of 2 appearing in front of the reaction term  aMFa 2 in Eq. (2) 4 accounts for the conservation of mass and equals the stoichiometric coefficient  = 2 of fatty acids [Fa] in the corresponding chemical reaction, Eq. (1) 1 .Further, in accordance with the experimental observations presented in Hermans (2017), the diffusion coefficient  is assumed to vary between the diffusion coefficients of the two physical material phases, i.e., the paint without any metal soap and the fully crystallised metal soap.This variation is quantified by the local volume fraction of the crystalline metal soap  as where  p and  s are the diffusion coefficients for fatty acid diffusion in the paint and in the crystalline metal soap, with the subscripts ''p' ' and ''s'' referring to ''paint'' and ''soap'', respectively.Further, the exponent  ≥ 1 is a constant that determines the rate at which the diffusion coefficient decreases as a function of the volume fraction  of crystalline metal soap formed.The value of  is determined by the ratio between the current concentration [cMFa 2 ] and the maximum concentration [cMFa 2 ] max of crystalline metal soap, i.e., The maximum concentration of crystalline metal soap [cMFa 2 ] max is dictated by the initial concentration of metal ions [M]( = 0) = [M] 0 , which, in accordance with the reaction equations given by Eq. ( 1), drives the formation of amorphous and crystalline metal soaps.Based on the reaction equations, Eq. ( 1), and considering that metal ions and metal soaps do not diffuse, from the conservation of mass the maximum concentration of crystalline metal soap can be computed as Eqs.
(2), to be completed with appropriate initial and boundary conditions, describe the concentration development of the metal ions, amorphous and crystalline metal soaps and free saturated fatty acids in the painting.According to these equations, metal soap formation can in principle occur at any location in the domain where free saturated fatty acids are present, provided that  ≠ 1.As discussed above, the metal soap growth process is assumed to depart from a crystalline nucleus of a small, specified size, for which  = 1.The chemical reactions are allowed to occur only in a small transition region located in between the metal soap crystal and the adjacent paint material, termed the reaction zone.Within this reaction zone, the volume fraction of crystalline metal soap is characterised by 0 <  < 1, and the evolution of all chemical species is governed by the full system of (diffusion-)reaction equations, Eq. ( 2).The reaction zone is defined via a scanning algorithm that identifies the uncrystallised material points present within a small, prescribed radial distance from the crystallised material points located closest to the current boundary of the metal soap crystal, see Eumelen et al. (2019) for more background information on the scanning algorithm.In the remaining domain, where either  = 0 (paint) or  = 1 (crystalline metal soap), the reaction term  cMFa 2 is set to 0 and the diffusion coefficient is specified, respectively, as  p (paint) or  s (crystalline metal soap).

Mechanical model 2.2.1. Continuum model
The mechanical model applied in this work is based on the formulation presented in Eumelen et al. (2019), as reviewed below.In any material point of the domain, the total strain tensor  is defined as the sum of the elastic strain tensor  e and a growth strain tensor  g : The growth strain tensor  g describes the volumetric expansion of the metal soap during the crystallisation process (Noble and Boon, 2007;Keune and Boon, 2007;Boon et al., 2003) and is taken to be proportional to the volume fraction  of crystalline metal soap (Eumelen et al., 2019): where  is defined in accordance with Eq. ( 7),  g is the growth strain associated to the formation of crystalline metal soap, and  is the second-order identity tensor.The formation and growth of crystalline metal soap thus induces stresses and strains in the paint layer, thereby generating a coupling between the chemical and mechanical fields.For material points in the paint phase ( = 0), the growth strain is equal to zero, and the total strain reduces to the elastic strain: Conversely, in material points located in the fully crystallised metal soap phase ( = 1), the growth strain reaches its maximum value  g s =  g , whereby the total strain is expressed by Finally, for material points in the reaction zone, where the paint and the metal soap phases are jointly present (0 <  < 1), Eq. ( 8) specifies into: where the subscript rz refers to ''reaction zone''.Departing from the above description, a constitutive assumption needs to be made for the two material phases considered, i.e., the paint and the crystalline metal soap.In general, the mechanical properties of oil paints depend on the type of pigment particles, drying oils and pigment to binder volume ratio, as well as the age of the paint and the specific environmental conditions, see Mecklenburg et al. (2004, 2012), Fuster-López et al. (2016), Hagan (2017), Fuster-López et al. (2019).In the present work, the ''bulk behaviour'' of the paint and the crystalline metal soap, for simplicity, is considered to be isotropic linear elastic, in accordance with the expressions where 4  p and 4  s are the fourth-order elastic stiffness tensors of the paint and crystalline metal soap, respectively, and  p and  s are the corresponding stress tensors.Further, the symbol ''∶'' denotes a double tensorial contraction.In order to describe the transformation from the paint material to the crystalline metal soap material, a basic phase transformation model is constructed as follows: Consider a material point located in the interfacial reaction zone, in which both the paint and the metal soap phases are present (0 <  < 1).Here, the total strain is assumed to be the same for both phases,  rz =  s =  p , i.e., Voigt's assumption is adopted, as a result of which the effective Cauchy stress tensor  rz follows from the weighted volume average of the stresses in the individual constituents (Eumelen et al., 2019;Hille et al., 2009): With the above constitutive models for the paint, crystalline metal soap, and their mixture in the reaction zone, Eqs. ( 13) and ( 14), in the absence of body forces the material response in any point of the domain is governed by the mechanical equilibrium equation Eq. ( 15) is completed by the appropriate boundary conditions.

Interface damage model
The stress field generated in a paint layer as a result of metal soap formation may trigger the onset and propagation of cracks.The discrete cracking behaviour of the paint is simulated by surrounding the continuum elements in the finite element model with interface elements equipped with a mixed-mode interface damage model.This modelling strategy was originally proposed in Xu and Needleman (1994), and allows for the robust simulation of crack patterns at arbitrary locations and in arbitrary directions.It naturally includes the effects of crack bifurcation, crack branching and crack coalescence, as previously demonstrated for applications related to historical paints (Eumelen et al., 2019(Eumelen et al., , 2020(Eumelen et al., , 2021)), wood (Luimes et al., 2018;Scheperboer et al., 2019;Luimes and Suiker, 2021), polymers (Tijssens et al., 2000), fibrous composites (Cid Alfaro et al., 2010a,b;Geng and Suiker, 2019) and cementitious materials (Scheperboer et al., 2021).The interface elements are characterised by the interface damage model proposed in Cid Alfaro et al. (2009).The fracture process in a material point initiates when the effective interfacial traction attains the value of the ultimate fracture strength of the material  u , which corresponds to a relative crack face separation  0 , as shown in Fig. 4. The fracture process is defined by a linear softening branch, for which the tensile strength progressively decays as a function of a damage parameter  that can vary between 0 (no damage) and 1 (damage completion).The elastic interfacial stiffness consequently decreases as (1 − ).The evolution of the fracture process is monitored by a damage history variable .The fracture process is completed when the interfacial material point has lost its strength ( = 1), which occurs when the relative ultimate crack face separation  u is reached.The toughness  c is obtained as the area under the traction-separation law, with  c =  u  u ∕2.A damage loading function is defined that governs the loading and unloading conditions during fracture.The evolution of the fracture process is described by a rate-dependent kinetic law.The adopted model additionally accounts for the mode-mixity of the fracture process, by considering a combination of the crack face separations in the normal and tangential directions of the crack.This enables to distinguish between the contributions of mode I (tension) and mode II (shear) cracking.For more details on the interface damage model and its numerical implementation, the reader is referred to Cid Alfaro et al. (2009).
The presence of cracks in the domain has consequences for the chemo-diffusive response of the system, since, across the crack surfaces, the diffusion of the free saturated fatty acids is reduced.In order to account for this coupling effect, following Eumelen et al. (2019), Luimes and Suiker (2021), in the diffusion-reaction model the concentration of the free saturated fatty acids across a discrete crack is assumed to be discontinuous while the flux is taken as continuous.Denoting the concentration jump of free saturated fatty acids across a crack as [[Fa ]], the corresponding flux is  Fa =  Fa ⋅ , where  Fa is the mass flux vector of the fatty acid species Fa, and  is the unit vector normal to the crack faces.The flux-concentration relation is given by where  c is the crack-reduced diffusion coefficient of the paint material, which is assumed to evolve as a function of the damage parameter  as

Numerical solution procedure
The coupled chemo-mechanical analysis is performed with the aid of the commercial finite element package ABAQUS standard,1 employing tailored user-subroutines for incorporating (i) the scanning algorithm described in Section 2.1 that is used for identifying and updating the interfacial reaction zone, (ii) the phase transformation model applied in the reaction zone, Eq. ( 14), (iii) the interface damage model employed for modelling discrete cracking, as outlined in Section 2.2.2, and (iv) the crack-reduced diffusion effect following from Eqs. ( 16) and ( 17).The chemical and mechanical field variables are solved for each time increment using a staggered, incremental-iterative update procedure, with the couplings between the chemical and mechanical fields established via a temporal extrapolation.The time increment is chosen relatively small, such that the error introduced by the time discretisation has a negligible influence on the numerical result.Each time increment starts with solving the coupled set of (diffusion-)reaction equations of the chemical model, as given by Eq. ( 2).The reactions between the various chemical species are hereby considered in the actual (updated) interfacial reaction zone between the paint material and the crystalline metal soap, which is constructed by applying the scanning algorithm described in Section 2.1.The flux of the free saturated fatty acids across crack surfaces is determined via Eq.( 16), whereby the crack-reduced diffusion coefficient in the paint material is determined from Eq. ( 17) using the damage value  computed at G.J.A.M. Eumelen et al. the previous time increment.The solution of the coupled diffusionreaction, Eqs. ( 2), provides the concentration fields of all chemical species, i.e., the concentrations of the free saturated fatty acids [Fa], metal ions [M], amorphous metal soap [aMFa 2 ] and crystalline metal soap [cMFa 2 ].With this result, the volume fraction of crystalline metal soap  is updated via Eq.( 7).The values of  are next transferred to the mechanical analysis, where they are used for calculating the growth strain tensor  g in accordance with Eq. ( 9).Additionally, the stresses in the paint and crystalline metal soap are computed from Eq. ( 13), and the stresses inside the interfacial reaction zone are calculated as a function of  via Eq.( 14).With this result, the equilibrium Eqs. ( 15) governing the mechanical problem are solved in an iterative fashion, thereby accounting for the nucleation and propagation of discrete cracks in accordance with the interface damage model outlined in Section 2.2.2.The converged iterative procedure provides the updates of the overall displacement field and the damage parameter  in the integration points.The updated damage parameter  is subsequently transferred to the chemical analysis, after which the incremental-iterative solution procedure above is repeated for the next time increment.

Parameter calibration for the diffusion-reaction model
The material parameters for the chemical model described by Eqs.(2), ( 5) and ( 6), are calibrated from the experimental measurements reported in Hermans (2017), which are briefly reviewed in Section 3.1.The calibration procedure applied and the parameter values obtained are presented in Sections 3.2 and 3.3, respectively.

Review of the experiment
The experiments discussed in Hermans (2017) focus on measuring the time evolution of the concentration of crystalline metal soap forming from a lead-containing ionomer exposed to a fatty acid solution.An ionomer sample was prepared that consisted of a lead sorbate 2 mixed with cold-pressed linseed oil.Note that recent experimental works on similar model systems containing zinc sorbate have revealed that zinc sorbate samples exhibit mechanical properties, solvent diffusion behaviour and a molecular structure that are comparable to those of typical zinc white oil paint films (Hermans et al., 2019;Baij et al., 2018b).The employed lead sorbate model system can thus be considered as representative of real oil paint layers.
The prepared mixture of linseed oil and lead sorbate was next spread on a glass plate and cured overnight in an oven.The resulting ionomer film, which was characterised by an average thickness  s = 150 μm, 2 A salt from sorbic acid.
was next placed at the bottom of an attenuated total reflection-Fourier transform infrared (ATR-FTIR) measurement cell.The measurement cell was subsequently filled with a solution of fatty acid (i.e., palmitic acid) and acetone over a height  Fa that ranged between 10 mm and 20 mm, as schematically indicated in Fig. 5(a).In the designed experimental setup, a one-dimensional diffusion-reaction process is activated, whereby the fatty acid diffuses along the vertical direction in the fatty acid-acetone mixture.Furthermore, the fatty acid, together with acetone, diffuses in the ionomer sample below, whereby it reacts to form metal soap.The process of metal soap formation was monitored by measuring the concentration of crystalline metal soap formed at the bottom of the ionomer sample, using time-dependent ATR-FTIR measurements.

Reference model and calibration procedure
In order to identify the kinetic reaction constants and the diffusion coefficients associated to the diffusion and reaction of saturated fatty acids in the ionomer film, the experimental set-up indicated in Fig. 5(a) has been idealised to a one-dimensional model.In this model, which is shown in Fig. 5(b), the height of the saturated fatty acid solution and the thickness of the ionomer sample are represented by two adjacent segments of length  Fa and  s , respectively.The thickness of the ionomer sample equals  s = 150 μm (Hermans, 2017).The height of the palmitic acid-acetone solution has been taken equal to the lowest value considered in the experiments,  Fa = 10 mm.Note however, that a variation of this value in the range of 10 to 20 mm is expected to have a negligible effect on the experimental result, since the sample thickness  s is much smaller than the height  Fa of the saturated fatty acid solution, i.e.,  s ∕ Fa ≪ 1.This hypothesis has been confirmed by simulations not presented here.
In the experiment, the saturated fatty acids diffuse through the acetone solution, and further diffuse and react within the ionomer below, see Fig. 5(b).Accordingly, the one-dimensional diffusion process in the acetone solution can be described via the diffusion equation: where  a is the diffusion coefficient of the saturated fatty acid in acetone, whereby the subscript ''a'' refers to the ''hosting medium'' (i.e., acetone).Further, for the ionomer sample it is assumed that the diffusion-reaction processes leading to metal soap formation are governed by the one-dimensional formulation of the system of chemical equations given by Eq. ( 2), complemented with the definitions in Eqs. ( 5), ( 6) and ( 7).The boundary value problem described by the geometry presented in Fig. 5(b) and governed by Eqs.(18) (acetone solution) and (2) (ionomer sample) must be completed by initial and boundary conditions.In the measurement cell, the initial concentration of saturated fatty acid is uniform and equal to [Fa]( = 0) = [Fa] 0 = 55 mol/m 3 (Hermans, 2017).In addition, zero-flux boundary conditions  Fa = 0 are adopted at the external boundaries  = 0 and  =  Fa +  s , which reflect that no saturated fatty acids were added to the palmitic acid-acetone solution after the experiment had started, and that the saturated fatty acids cold not diffuse outside the measurement cell.Furthermore, the initial metal ion concentration [M]( = 0) = [M] 0 is taken as uniform within the ionomer sample and equals [M] 0 = 160 mol/m 3 (Hermans, 2017). Eqs.
(2), ( 5), ( 6), ( 7) and ( 18) have been implemented in the mathematical software program MATLAB, and solved using the standard partial differential equation solver.The computational result has been subjected to a non-linear, least-squares fitting procedure, in order to match the concentration of crystalline metal soap [cMFa 2 ] measured at the bottom of the sample,  =  Fa +  s , as reported in Hermans (2017).For this purpose, the diffusion coefficient of the saturated fatty acid in acetone has been selected as  a = 6.1 ⋅ 10 −11 m 2 /s (Hermans, 2017).Correspondingly, the calibrated parameters are the reaction rates  + ,  − and  crys , the diffusion coefficients  p and  s , and the exponent  in Eq. ( 6) that governs the transition of the diffusion coefficient in the ionomer from  p to  s .

Result of calibration procedure
Fig. 6 shows the time evolution of the concentration of crystalline metal soap [cMFa 2 ], evaluated at the bottom of the ionomer sample, as obtained from the experimental measurement (black line) and the numerical simulation (red line).From the solution of Eqs. ( 2) and ( 18), combined with the calibration procedure described above, the model parameters have been determined in correspondence with an  2 -value of 0.9748.Accordingly, the kinetic reaction constants are obtained as  + = 0.491⋅10 −6 m 6 ⋅ (mol 2 s) −1 ,  − = 0.169 s −1 and  crys = 0.600 s −1 .The diffusion coefficients in the paint and in the metal soap are  p = 5.100 ⋅ 10 −11 m 2 /s and  s = 3.635 ⋅ 10 −16 m 2 /s, respectively, and the exponent defining the transition in the ionomer diffusion coefficient is  = 2.The calibrated parameters are summarised in Table 1.It is noted that the above calibration procedure somewhat differs from the calibration procedure applied in Hermans (2017).For most of the chemical model parameters this difference in the calibration procedure had a minor effect on the actual values obtained, as it should, with the exception being the value of  s in Eq. ( 6).This difference, however, is consistent with a somewhat lower exponential value of  = 2 in Eq. ( 6).It may thus  2) for the ionomer sample and Eq. ( 18) for the acetone-fatty acid solution.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)be concluded that the calibrated parameters realistically represent the diffusion-reaction behaviour in the experimental system.Nevertheless, it should be mentioned that the presence of acetone in the experimental system promotes diffusion of saturated fatty acid in the ionomer film as well as swelling of the layer (Baij et al., 2018a), such that the calibrated diffusion parameters of the fatty acid in the ionomer overestimate those in a real historical painting.In order to quantitatively improve the model predictions, future research should aim at refining the estimated values of the input parameters by dedicated non-invasive experiments, or, when accurate experimental values are impossible to obtain, at investigating the sensitivity of the overall system behaviour to key parameter values or mechanistic assumptions.

Numerical example 1: Single-layer paint system
In this section a numerical example is analysed that is representative of degradation phenomena in oil paintings associated to the formation of multiple metal soap protrusions, as observed in historical paintings,

Table 1
Diffusion coefficients and reaction rates calibrated from the experiment reported in Hermans (2017) in which a palmitic acid diffuses and reacts in a lead-containing ionomer, in accordance with the representation in Fig. 5 2019) that neglects the contribution of the amorphous metal soap in the chemical reaction process.Hence, a comparison between the results obtained with the original model and the current, enhanced model described in Section 2.1 allows to reveal the important contribution of the amorphous metal soap on the degradation process.Section 4.1 describes the geometry of the paint layer model, the material properties, and the initial and boundary conditions.The results of the numerical simulation are presented in Section 4.2.

Geometry and FEM discretisation
The single-layer paint system is modelled as a two-dimensional rectangular domain with a thickness  = 15 μm and a width 4 = 60 μm, in which three metal soap nuclei are present, see Fig. 7.The shape of the nuclei is taken as circular, with an initial radius equal to one tenth of the layer thickness,  0 = ∕10 = 1.5 μm.The finite element discretisation of the geometry is similar to that used in Eumelen et al. (2021), whereby, after a mesh convergence study, the element size close to the metal soap nuclei was set to 0.2 μm, and at the corners of the geometry to 0.5 μm, as required for obtaining mesh-independent numerical results.Accordingly, the computational domain is discretised with 27, 030 continuum elements.For the diffusion-reaction analysis, 3node isoparametric elements with a 3-point Gauss quadrature are used, while for the mechanical analysis 6-node plane-strain isoparametric elements with a 3-point Gauss quadrature are employed.The interface elements describing the discrete fracture processes in the FEM geometry are included by using a dedicated Python script, which first reduces the size of all the original continuum elements by 1%, then detects the opposing element edges, and subsequently places interface elements in between the continuum elements, leading to 40, 345 interface elements in total.For the chemical analysis 4-node interface elements with a 2-point Gauss quadrature are used, and for the mechanical analysis 6-node interface elements with a 3-point Lobatto quadrature are selected.Note that the element formulation for the mechanical interface elements is not available in ABAQUS Standard, and is therefore incorporated by means of a user-defined element routine based to the numerical implementation presented in Schellekens and de Borst (1993).Finally, the radial detection distance used in the scanning algorithm for defining the thickness of the interfacial reaction zone from the integration points closest to the boundary of the crystalline metal soap is selected to be 0.2 μm, which corresponds to the smallest characteristic element size in the FEM model.Hence, the radial detection distance is 7.5 times smaller than the initial radius  0 of the metal soap nucleus, and 75 times smaller than the thickness  of the paint, and has a negligible influence on the computational results.

Material properties
The paint layer is assumed to be made of lead white paint.For the diffusion-reaction analysis the material properties obtained from the calibration procedure discussed in Section 3 have been used, as summarised in Table 1.The input parameters for the mechanical analysis are presented in Table 2.The Young's modulus of the lead white paint is extracted from the uniaxial stress-strain data presented in Fuster-López et al. ( 2016), leading to  p = 115 N/mm 2 , and the Poisson's ratio is set to  p = 0.3 (Eumelen et al., 2019).For the crystalline metal soap, no experimental data on the mechanical properties are available in the literature.Nevertheless, in Eumelen et al. (2019) it was demonstrated that the value of the mismatch between the stiffness of the metal soap crystal and that of the oil binder influences the growth characteristics of the metal soap crystal only marginally.Hence, for simplicity the elastic constants for the metal soap crystal are taken equal to those of the paint material, i.e.,  p =  s = 115 N/mm 2 and  p =  s = 0.3.The fracture properties of the paint and crystallised metal soap are derived from the stress-strain diagram reported in Fuster-López et al. ( 2016), leading to an ultimate tensile strength  u 1 = 2.2 N/mm 2 and a mode I fracture toughness  I,c = 0.418 N/mm.From the simulation results presented in Eumelen et al. (2021), it can be observed that the fracture process in the painting is indeed dominated by mode I cracking events.Therefore, the fracture properties in mode II are expected to only have a minor influence on the computational results, and, for simplicity, are taken equal to the properties under mode I conditions, i.e.,  u 2 =  u 1 = 2.2 N/mm 2 and  II,c =  I,c = 0.418 N/mm.The elastic stiffness of the interface elements is set to a relatively high value,  = 10 6 N/mm 3 , which, together with the choice of an appropriate mesh fineness, ensures that artificial response contributions related to the use of an elastic interface stiffness  in the traction-separation law are negligible, see Cid Alfaro et al. (2010a) for more details.Accordingly, the elastic response of the layer model is correctly determined by the continuum elements through the elastic properties of the paint and the crystalline metal soap.
To the best of the authors' knowledge, for the chemical growth strain associated to metal soap formation no experimental value has been reported in the literature.Hence, the value of  g * = 0.1 pragmatically adopted for the plane-stress analyses presented in Eumelen et al. (2019Eumelen et al. ( , 2020Eumelen et al. ( , 2021) ) is adjusted for the current, plane-strain analysis by requiring that the effective, in-plane chemical growth strain in both 2D models is the same, which, in accordance with Equation (26) in Eumelen et al. (2019), is the case for  g = (1 +   ) g * = (1 + 0.3) ⋅ 0.1 = 0.077.Note hereby that the choice of representing a 3D paint layer by a 2D plane-stress or plane-strain model is somewhat arbitrary, but, in principle, leads to comparable results for the in-plane chemo-mechanical behaviour of the paint layer.

Initial and boundary conditions
In the chemical model the initial concentration of saturated fatty acid is assumed to be uniform throughout the domain.The saturated fatty acid concentration has been calculated based on the composition of the lead white paint analysed in Fuster-López et al. ( 2016), which contains stearic acids3 and lead-based pigments.The initial fatty acid concentration [Fa]( = 0) = [Fa] 0 is taken as maximal by assuming that all stearic acids from the binder are available for the formation of metal soap, in accordance with dividing the volumetric weight of the fatty acid in the oil binder by its molecular weight.The volumetric weight of the fatty acid is calculated by multiplying the volume fraction of linseed oil in a lead white paint (=0.556 (Fuster-López et al., 2016)) by its volumetric weight (=930 kg/m 3 (Alfred, 2000)) and by the weight fraction of stearic acid in linseed oil (=0.05 (Alfred, 2000)).By adopting a molecular weight for the stearic acid of 0.2835 kg/mol (Anneken et al., 2012), an initial fatty acid concentration of [Fa] 0 = 91 mol/m 3 is obtained.
The initial concentration of lead ions [M]( = 0) = [M] 0 is determined by assuming that only a specific fraction  of the lead ions present in the pigment particles reacts to form crystalline lead soap, see also the discussion below Fig. 3 in Section 1.For this purpose, consider the density of 6700 kg/m 3 of the lead white pigment (Lide, 2003) and the density of 1403 kg/m 3 of the lead stearate (i.e., lead soap) (Martínez-Casado et al., 2017).The molecular densities of lead white pigment and lead stearate follow from dividing these densities by the corresponding molecular weights, i.e.,  pigm = 6700∕0.776= 8634 mol/m 3 for the lead white pigment and  s = 1403∕0.774= 1812 mol/m 3 for the lead stearate.From the chemical formula of lead white pigment, 2PbCO 3 ⋅Pb(OH) 2 , it can be concluded that one mole of lead white pigment contains three moles of lead ions.Hence, the total molecular density of lead ions in lead white pigment is  Pb,pigm = 25902 mol/m 3 .Consider now a unit volume of original lead white paint material  paint,0 = 1 m 3 , which contains a volume fraction of 0.444 of lead white pigment (Fuster-López et al., 2016) and a volume fraction of 1 − 0.444 = 0.556 of binder.The total number of moles of lead ions in the lead white pigment then becomes  Pb,pigm = 0.444 × 25902 = 11500 mol.After the formation of lead soap, the new paint volume  paint is obtained as the sum of the volume of remaining lead ions in the lead white pigment  Pb,pigm , the volume of lead soap  s , and the volume of the original binder material  bin (=0.556 m 3 ): The volume occupied by the lead soap follows from where  is the fraction of lead ions available for the reaction into lead soap.The volume of the remaining lead ions in the paint then equals The new paint volume can be calculated from the volumetric expansion induced by the free chemical growth strain  g * as which, with  g * = 0.1 (Eumelen et al., 2019(Eumelen et al., , 2020(Eumelen et al., , 2021) ) and an original paint volume of  paint,0 = 1 m 3 , leads to  paint = 1.3 m 3 .Combining this result with Eqs. ( 19) to ( 21), and solving for , results in  = 0.051.Hence, the initial concentration of lead ions in 1 m 3 lead white paint that is available for the reaction into lead soap becomes In the simulations, it is realistic to assume that the initial concentration of lead ions is zero within the crystalline lead soap nuclei.In addition, the initial concentration of amorphous lead soap in the paint material is taken as zero, [aMFa 2 ]( = 0) = [aMFa 2 ] 0 = 0 mol/m 3 .
The chemical boundary conditions required for the diffusionreaction analysis need to be prescribed in terms of the diffusing species, i.e., the saturated fatty acids.Accordingly, at the left and right domain boundaries the value of the initial concentration [ Fa] = [Fa] 0 is imposed, see Fig. 7. Furthermore, at the top and bottom boundaries zero-flux boundary conditions are applied,  Fa = 0.At the bottom of the paint layer this boundary condition reflects that the free saturated fatty acids do not penetrate the substrate below.
In the mechanical model the paint layer is initially considered to be stress-free.At the bottom boundary the vertical displacement is constrained, which represents the relatively stiff support provided to the paint by the wooden or canvas substrate.At the left and right boundary of the domain the horizontal displacement is prescribed to be zero.This boundary condition represents the upper limit for the horizontal deformation constraint provided by the paint material outside the computational domain.The top boundary of the domain is stress-free, in correspondence with zero-traction boundary conditions.

Results
For generality, the computational results are presented in terms of dimensionless parameters, with the dimensionless time defined as  =  p ∕ 2 , where  p and  are the diffusion coefficient in the paint material and the thickness of the paint layer, respectively.Since cracking in the paint occurs predominantly under mode I conditions, the local amount of cracking is assessed in terms of the normal (mode I) crack opening  1 , which in dimensionless form reads v1 =  1 ∕.thereby showing the progressive growth of the three metal soap nuclei and the formation of discrete cracks.The contour plot variable indicates the volume fraction  of crystalline metal soap, with the red colour ( = 1) referring to fully crystallised metal soap and the blue colour ( = 0) designating the original paint material.The narrow reaction zones defining the transition from the paint material to a metal soap crystal are shown by the range of colours related to 0 <  < 1.The growth of the metal soap crystals and the associated volumetric expansion clearly cause an upward deflection of the upper surface of the paint layer.This is particularly evident at the left crystal, which is able to reach the upper part of the paint layer and then erupts through the outer surface.It can be further seen that the three crystals initially grow at approximately the same rate and develop a circular shape, see Figs. 8(a), (b) and (c).Note that this shape closely resembles that of the metal soaps visible in the cross-sectional image taken from the ceiling painting in the Room of Trustees, Burgerweeshuis, Zierikzee, see Fig. 1(a).Nevertheless, at the largest time instant considered, t = 10185, the central and right crystals have coalesced into a single, elongated metal soap aggregate, see Fig. 8(d).This elongated shape shows resemblance with that of the metal soap crystals observed in View of Delft by Johannes Vermeer, see Fig. 1(b), although the approximate crystal size of 100 μm observed in this painting is somewhat larger than the maximal crystal size of 25 μm following from the numerical simulation.Trivially, the maximum crystal size is determined by the specific time frame and the paint layer geometry considered in the analysis.
During the process of metal soap growth, discrete cracks are initiated once the stress in the paint material attains the fracture strength.As illustrated in Fig. 8(c), near the top surface of the paint layer, just above the left metal soap crystal, 6 distinct mode I cracks develop.
The time evolution of the dimensionless normal crack opening v1 =  1 ∕ of the 4 largest mode I cracks is shown in Fig. 9(a), and is evaluated at the upper surface of the paint layer at which the crack mouth opening is maximal.As shown in Figs.8(b) and (c), the first surface crack nucleates close to the vertical symmetry line of the left metal soap crystal, and partly penetrates the crystal.The additional surface cracks develop in the paint material region adjacent to the first surface crack, see Fig. 8(c), thereby inducing a further relaxation of the tensile stresses at the paint surface.When the left metal soap crystal approaches the tips of these surface cracks, the cracks start to close, due to compressive stresses generated in the crystal by restrained volumetric expansion originating from the chemical growth strain, see Eq. ( 9).As illustrated in Fig. 9(a), the completion of crack closure occurs at approximately the same time instant for all 4 surface cracks.Note further from Figs. 8(c) and (d) that similar mode I surface cracks develop above the right metal soap crystal, although later in time.In addition, diffusive patterns of micro-cracks are generated in between the three metal soap crystals, which occurs along trajectories directing towards their geometrical centres, see Figs. 8(b), (c) and (d).As for the surface cracks, crack closure is observed once the metal soap crystals approach these micro-cracks, see Fig. 8(d).
In order to reveal the specific influence of the enhanced chemical model on the growth processes of the crystalline metal soaps and the discrete cracks, the present numerical results are compared to those in Eumelen et al. (2021) that are based on the original model proposed in Eumelen et al. (2019), which, as mentioned, neglects the reversible reaction of amorphous metals soaps and only accounts for the formation of irreversible, crystalline metal soap.The numerical simulation presented in Eumelen et al. ( 2021) starts from the same Fig. 9. Single-layer paint system: Time evolution (using the dimensionless time t) of the normalised crack opening  1 of 4 dominant mode I surface cracks that develop above the left metal soap crystal, as predicted by (a) the current, enhanced model that includes the formation of amorphous metal soap, and (b) the original model presented in Eumelen et al. (2019Eumelen et al. ( , 2021) ) that neglects the formation of amorphous metal soap.The vertical dashed line in Fig. 9(a) designates the specific time t = 5500 at which the left metal soap crystal shown in Fig. 8 reaches a relative size of r = ∕ = 0.27, which corresponds to the final crystal size reached at t = 6.5 in the simulation considered in Fig. 9(b).initial configuration as that displayed in Fig. 8(a), and shows that the growth rate and final size and shape of the three metal soap crystals strongly depend on their precise location in the computational domain.Specifically, the two outer metal soap crystals eventually obtain a slightly oval-shaped geometry that is almost two times larger in size than that of the central, circular metal soap crystal, see Figure 4 in Eumelen et al. (2021).This result is rather different from that shown in Fig. 8, and can be ascribed to the fact that, by ignoring the intermediate amorphous metal soap phase, in the original model the free saturated fatty acid is consumed at a significantly higher rate, whereby its spatial distribution, and thus the corresponding development of metal soaps crystals, becomes non-uniform across the paint layer, see Figure 6 in Eumelen et al. (2021).Conversely, the relatively low rate at which the free fatty acid is consumed in the current, enhanced model ensures that its spatial distribution remains uniform when time develops, thus resulting in a virtually homogeneous growth of the three individual metal soap crystals (as long as they do not coalesce).
The fracture response in Fig. 9(a) is now compared to that in Fig. 9(b), which illustrates the time evolution of the dimensionless normal crack opening v1 =  1 ∕ of the 4 main surface cracks as computed by the original model that neglects amorphous metal soap formation.The 4 surface cracks considered in Fig. 9(b) are located above the left metal soap crystal, at comparable locations as shown in Fig. 8 for the enhanced model, and appear to have similar crack widths as those in Fig. 9(a).The simulation carried out with the original model covers a total time of t = 6.5, at which the normalised average radius of the left metal soap crystal has reached a value of  = ∕ = 0.27.For a proper comparison with the enhanced model, in Fig. 9(a) the simulation time corresponding to this crystal size has been indicated by a vertical dashed line.It can be observed that the vertical dashed line is related to a dimensionless time of  = 5500, which is about a factor of 850 larger than the dimensionless time t = 6.5 at which the crystal size  = ∕ = 0.27 is reached in the original model.In addition, Fig. 9(a) illustrates that for the enhanced model the time at which the first surface crack nucleates is approximately 1300 times larger than for the original model in Fig. 9(b).This again confirms that the addition of the amorphous metal soap phase in the enhanced model substantially delays the formation of crystalline metal soap, and thus the generation of discrete cracks, and shows that the characteristic time scale associated to the chemo-mechanical degradation of the single paint layer is increased by approximately three orders of magnitude.Specifically, from the definition of the dimensionless time  =  p ∕ 2 , it follows that in the enhanced model the metal soap formation process occurs over a time frame of hours to days, while in the original model it takes place over a time frame of seconds to minutes.
Although the characteristic time scale computed by the current, enhanced model is in agreement with that of the experiments reported in Hermans (2017), see also Fig. 6, it is still significantly smaller than the time scale of years to decades characterising metal soap degradation processes of historical paintings in museum collections.This difference can be ascribed to the relatively high concentration value used in the definition of the initial and boundary conditions for the saturated fatty acid.As described in Section 4.1.3,the initial and boundary conditions imposed are representative of the maximum saturated fatty acid concentration available for the chemical reaction process; in specific, it is assumed that all stearic acids from the binder are available for the formation of metal soap.It is further assumed that this maximum concentration is already present in the paint layer from the beginning of the simulation.In real, historical paintings, however, the concentrations of saturated fatty acids and metal ions gradually build up in paint layers as a result of the drying process of the paint, which typically takes place over a period of years to decades (Baij et al., 2019;Hermans et al., 2015Hermans et al., , 2016a;;Baij et al., 2018a;Mecklenburg et al., 2004).The analysis of the dependency of the characteristic time scale of the chemo-mechanical degradation process on the drying process of paint is considered to be a topic for future research.

Numerical example 2: Multi-layer paint system
In this section the influence of the formation of a single lead soap crystal on cracking and delamination in a multi-layer paint system is analysed.The geometry, material properties, and initial and boundary conditions of the multi-layer model are presented in Section 5.1, and the results of the numerical simulations are discussed in Section 5.2.

Geometry and FEM discretisation
The geometry of the multi-layer paint system is characterised by a 2 × 2 square domain whereby  = 50 μm, which is composed of a ground layer indicated as layer 1, and a pictorial layer indicated as layer 2, see Fig. 10.The ground layer is made of lead white paint and has a thickness  1 = 85 μm, while the pictorial layer consists of cobalt blue paint and has a thickness  2 = 15 μm.Note that the cobalt blue paint has been chosen as a representative paint material that does not contain Fig. 10.Geometry and essential chemical and mechanical boundary conditions for the multi-layer paint system.In the ground layer (layer 1, grey colour) lead ions are present, while in the pictorial layer (layer 2, white colour) no lead ions are present.The red circle indicates the lead soap nucleus initially present in the ground layer.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)lead ions.The ground layer contains a single, circular lead soap nucleus with an initial radius  0 = 1.5 μm, which is located along the vertical symmetry axis of the domain, at a distance  n = 60 μm from the bottom of the paint system.The modelled domain is discretised by using 44,856 continuum elements and 67,185 interface elements.The element types applied in the simulation are the same as those selected for the singlelayer model, see Section 4.1.1 for more details.The minimum element size is applied close to the lead soap nucleus and equals 0.2 μm.The maximum element size taken at the left and right domain boundaries is 5 μm.The radial detection distance used in the scanning algorithm for defining the thickness of the interfacial reaction zone between the lead soap crystal and the paint material is 0.4 μm.It has been confirmed that for the chosen mesh density the numerical results converge within a small, specified tolerance.

Material properties
The material properties in the chemical model are presented in Table 1, as obtained from the calibration procedure presented in Section 3.For simplicity, the diffusion coefficients of the free saturated fatty acids in the ground layer and in the pictorial layer are taken the same.
The material parameters used in the mechanical model are listed in Table 3.For the ground layer (layer 1), the elastic modulus is derived from the uniaxial stress-strain response of an 18.75 year old lead white paint (Mecklenburg et al., 2012), leading to  p = 130 N/mm 2 .For the pictorial layer (layer 2), the elastic modulus is obtained from the stress-strain response of an 18 year old cobalt blue paint presented in Fuster-López et al. (2019), resulting in  p = 35 N/mm 2 .The Poisson's ratio is taken the same for the two paint layers and equals  p = 0.3.Since it is assumed that there are no lead ions present in layer 2, lead soap growth can only occur in layer 1, whereby the elastic properties of the lead soap crystal are taken the same as those of the lead white paint,  s =  p = 130 N/mm 2 and  s =  p = 0.3.The ultimate tensile strength of the lead white paint is obtained from the uniaxial stressstrain data presented in Mecklenburg et al. (2012), resulting in  u 1 = 4.5 N/mm 2 .The mode I toughness of the lead white paint is derived from the area under the uniaxial stress-displacement response -with the displacement determined as the strain multiplied by the gauge length  g = 76 mm (M.Mecklenburg, personal communication) -, leading to  I,c = 7.86 N/mm.Similarly, the ultimate tensile strength and mode I  2019), leading to  u 1 = 0.35 N/mm 2 and  I,c = 0.033 N/mm, with the latter value calculated using the applied gauge length of  g = 7.5 mm.For simplicity, for both layers the mode I and mode II fracture properties are taken the same,  u 2 =  u 1 and  II,c =  I,c .As for the single-layer paint model analysed in Section 4, the growth strain is selected as  g = 0.077 and the elastic stiffness in the interface damage model equals  = 10 6 N/mm 3 .Finally, the fracture properties of the delaminating cracks that may develop along the material interface between the ground layer and the pictorial layer are assumed to be equal to the fracture properties of the pictorial layer (layer 2), which is the weaker, more brittle layer.

Initial and boundary conditions
The initial conditions and boundary conditions adopted for the multi-layer paint model are comparable to those applied for the singlelayer paint model discussed in Section 4.1.3.In the chemical model the initial concentration of free saturated fatty acid is assumed to be uniform throughout the domain and equals [Fa] 0 = 91 mol/m 3 .The initial concentration of lead ions is taken as uniform in the lead white ground layer (layer 1), and equals [M] 0,1 = 587 mol/m 3 .As mentioned, no lead ions are present in the cobalt blue pictorial layer (layer 2), in correspondence with [M] 0,2 = 0.The initial concentrations of the amorphous and crystalline metal soaps are assumed to be equal to zero, i.e., [aMFa 2 ] 0 = [cMFa 2 ] 0 = 0 mol/m 3 .The above value for the initial concentration of free saturated fatty acid [ Fa] = [Fa] 0 is also imposed at the left and right boundaries of the computational domain.Furthermore, at the top and bottom boundaries of the domain zero-flux boundary conditions are adopted for the saturated fatty acid,  Fa = 0.
Similar to the single-layer paint model studied in Section 4, in the mechanical model the multi-layer paint system is initially considered to be stress-free.The vertical displacement is prevented at the bottom boundary of the domain, while at the left and right boundaries the horizontal displacement is constrained.Finally, the top boundary is stress-free, in correspondence with zero-traction boundary conditions.

Results
Fig. 11 shows the progressive growth of the lead soap crystal and the cracks and interfacial delamination generated in the multi-layer paint system, considering 4 different values of the dimensionless time  (with the deformation plotted at true scale).The dimensionless time is defined as  =  p ∕ 2 , with  = 50 μm taken as half of the total model thickness of the paint system, see also Fig. 10.The contour plot variable indicates the volume fraction of crystalline lead soap , with the red colour referring to the fully crystallised lead soap ( = 1) and the blue colour reflecting the original paint material ( = 0).The top grey layer indicates layer 2, in which no lead ions are present, and therefore does not contribute to the formation of lead soap.In addition, the intact interface between the ground layer (layer 1) and the pictorial layer (layer 2) is designated by the black horizontal lines, and the white lines in Figs.11(c) and (d) indicate two symmetrical, delaminating cracks emerging at the layer interface.
The lead soap crystal develops with a circular shape in the ground layer, see Figs. 11(a) and (b), and generates a mode I surface crack in the upper, pictorial layer when approaching the interface between the two layers, see Figs. 11(c) and (d).The volumetric expansion caused by the chemical growth strain induces an upward deflection of the upper, free surface of the paint system.In combination with the mismatch in the material properties of the ground layer and pictorial layer, relatively large shear stresses develop along the layer interface, which result in the nucleation and growth of two symmetrical, mode II delaminating cracks (indicated by the white lines).It can be seen from Fig. 11(d) that the delaminating cracks are close to coalescence once the mode I surface crack has grown through the entire thickness of the upper, pictorial layer, thereby forming a T-shape crack with the mode I surface crack.T-shape cracks are regularly observed in historical paintings, and may originate from metal soap formation and/or indoor climate fluctuations, often leading to local paint flaking and paint loss (Bosco et al., 2020(Bosco et al., , 2021)).
Fig. 12(a) presents the time evolution of the normalised normal crack opening v1 =  1 ∕ of the mode I crack surface crack, as evaluated at the surface of the upper, pictorial layer at which the crack mouth opening is maximal.It can be seen that the crack grows at an approximately constant rate.Contrary to what occurs in the single-layer paint system, see Fig. 9(a), crack closure does not occur within the time frame considered.This is due to the fact that the compressive stress caused by the constrained volumetric expansion of the lead soap crystal has partially relaxed as a result of the formation of the two mode II delaminating cracks at the layer interface, and that the lead soap crystal is not able to penetrate the upper pictorial layer due to the absence of lead ions in this layer.Fig. 12(b) illustrates the time evolution of the two mode II delaminating cracks that develop along the layer interface, expressed in terms of the normalised delamination length,  d =  d ∕ with  d the actual delamination length.Observe that the delaminating cracks grow in a synchronous (symmetrical) fashion at an initially high rate, but that the growth rate decreases substantially when time develops.As illustrated in Figs.12(c) and (d), the decrease in growth rate may be caused by the fact that the delaminations get closer to the mode I surface crack with increasing time, thereby reaching a region in which the stress cannot grow, such that the driving force for further delamination decreases.
It is finally noted that the small, staircase increments characterising the fracture responses in Fig. 12 result from the discrete interplay between the local element size and the radial detection distance used in the scanning algorithm for defining the thickness of the interfacial reaction zone.Indeed, in the limit of an infinitesimal element size and an infinitesimal reaction zone, the fracture responses in Fig. 12 are expected to become ideally smooth.

Conclusions
This paper proposes a chemo-mechanical model for the prediction of degradation phenomena in historical oil paintings as induced by metal soap reactions, such as the formation of metal soap aggregates, the deformation of the paint surface, paint cracking, and delamination between paint layers.Departing from the framework presented in an earlier work (Eumelen et al., 2019), the chemical model has been substantially enhanced to describe both the reversible reaction between free saturated fatty acids and metal ions forming amorphous metal soap, and the irreversible reaction into crystalline metal soap.The chemical parameters of the enhanced chemo-mechanical model have been calibrated on the experimental results published in Hermans (2017) that measure the time evolution of the concentration of crystalline metal soap forming from a lead-containing ionomer exposed to a fatty acid solution.Metal soap crystallisation introduces stresses in the paint that promote crack nucleation and propagation; this discrete cracking behaviour of the paint is simulated by surrounding the continuum elements in the finite element model with interface elements equipped with a mixed-mode interface damage model, which allows for the robust simulation of crack patterns at arbitrary locations and in arbitrary directions.
The enhanced chemo-mechanical damage model has been applied for analysing two types of boundary value problems, which capture the essential degradation mechanisms observed in historical oil paintings in museum collections.The first boundary value problem concerns the deformation and fracture response of a single-layer paint system under the growth of multiple metal soap crystals.The computational results show that the growth of the metal soap crystals and the associated volumetric expansion cause a substantial upward deflection of the top surface of the paint layer and induce multiple cracks at the layer surface and in between the individual metal soap crystals.The growth rate of the three metal soap crystals is independent of their location within the domain, which is due to the homogeneous spatial development of the fatty acid concentration in time.Further, at some stage two different metal soap crystals coalesce into a single, oval-shape metal soap aggregate that is characteristic for metal soaps observed in real, historical paintings.
The modelling results are compared to those presented in Eumelen et al. (2021) that were computed by the original chemo-mechanical model presented in Eumelen et al. (2019), which neglects the intermediate, amorphous metal soap phase and only accounts for the irreversible chemical reaction into crystalline metal soap.From this comparison, it follows that the inclusion of the intermediate, amorphous metal soap phase in the model formulation increases the characteristic time scale associated to the chemo-mechanical degradation process by three orders of magnitude, which is necessary to realistically predict the time-dependency of degradation processes observed in historical paintings in museum collections.Also, local differences are observed in the metal soap growth characteristics and crack profiles predicted by the two models, which are caused by the fact that with the current, enhanced chemo-mechanical model the saturated fatty acid profile develops homogeneously across the paint layer, while with the original chemo-mechanical model this occurs in a heterogeneous fashion.
The second boundary value problem considers the deformation, cracking and delamination processes in a multi-layer paint system, as induced by the growth of a single lead soap crystal.The multi-layer paint system is composed of a pictorial layer that is supported by a ground layer containing a lead soap nucleus.The growth of the lead soap nucleus induces a significant upward deformation of the paint surface.In addition, a mode I surface crack develops across the thickness of the upper, pictorial layer, which coalesces with two symmetrical delaminations developing along the layer interface, thereby forming a T-shape crack that ultimately may lead to local paint flaking and paint loss.
The enhanced chemo-mechanical formulation and modelling results presented in this communication provide a better understanding of the coupled chemo-mechanical degradation behaviour of historical oil paintings due to metal soap formation.Future research will aim at refining the model input parameters by, for example, accurately determining the mechanical parameters of small paint samples taken from historical paintings via nanoindentation testing (Salvant et al., 2011;Freeman et al., 2019;Fujisawa and Łukomski, 2019;Tiennot et al., 2020;Eumelen et al., 2022).Additionally, the effects of conservation treatments (Baij et al., 2019), indoor climate conditions (Keune et al., 2016), and paint drying characteristics (Mecklenburg et al., 2004) on the process of metal soap degradation need to be explored in more detail.The interaction between metal soap-induced degradation mechanisms and climate-induced failure of historical paintings (Bosco et al., 2020(Bosco et al., , 2021) ) is another topic for future study.Such investigations eventually should lead to optimal preventive conservation measures and environmental display/storage conditions required for preserving valuable historical paintings for future generations.

Fig. 1 .
Fig. 1.Examples of metal soap-induced degradation mechanisms.(a) Ceiling painting with allegorical representations, 1680, Room of Trustees, Burgerweeshuis, Zierikzee, The Netherlands (Photo: Wim Ruigoord).The detail shows a dark-field image of a sample cross-section taken from the red sky of Justitia and Mars in a Warm Embrace, right painting.In the bottom orange-red paint layer, which contains lead white and red lead, three lead soap protrusions have formed.An interface undulation can be observed in this cross-section, whereby the crystals in the bottom layer start to penetrate the layers above(Thoury et al., 2019).(b) Johannes Vermeer, View of Delft, ca.1660-1661, Oil on canvas, Mauritshuis, The Hague, The Netherlands.The detail shows a bright-field image of a paint cross-section taken at the left edge of the painting, where the green (now blue) bush is painted over the red-tiled roof, showing the formation of large metal soap protrusions in the dark-red paint layer that have erupted through the blue-green surface paint layer (van Loon, 2009; van Loon et al., 2012).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Schematic representation of the processes governing metal soap formation.(a) During the drying process of the paint film, the oil binder polymerises into a densely cross-liked polymer network by means of autoxidation reactions and releases free saturated fatty acids as a degradation product.Metal ions, which are released from the pigment particles, become complexed with the polymer network, thereby forming an ionomer.(b) Metal ions from the ionomer react by means of a reversible reaction with the free saturated fatty acid molecules to form amorphous metal soap.(c) Amorphous metal soap crystallises by means of an irreversible reaction and grows in large metal soap aggregates.

Fig. 3 .
Fig. 3. Transmission electron microscope (TEM) bright field image of a metal soap nucleus, which contains crystalline zinc soap dendrites (indicated in red) that are surrounded by paint material (indicated in blue) based on a linseed oil polymer.The picture (without the colours) is taken from Hermans et al. (2018).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. (a) Schematic representation of the experiment reported inHermans (2017).The experimental set-up consisted of a lead-based ionomer sample with a thickness  s located at the bottom of an ATR-FTIR measurement cell, which was filled with a mixture of palmitic acid and acetone over a height  Fa .The concentration of crystalline metal soap forming from the reaction of the palmitic acid with the lead ions from the ionomer was evaluated at the bottom of the ionomer by means of ATR-FTIR measurements.(b) One-dimensional schematisation of the experimental setup presented in Hermans (2017), with the corresponding initial and boundary conditions.

Fig. 6 .
Fig.6.Parameter calibration for the diffusion-reaction model.Time evolution (in minutes) of the concentration of crystalline lead soap (mol/m 3 ), as evaluated at the bottom of an ionomer sample,  =  Fa +  s .The black line represents the experimental result reported inHermans (2017), and the red line reflects the numerical result for the one-dimensional model in Fig.5(b) by solving Eq. (2) for the ionomer sample and Eq.(18) for the acetone-fatty acid solution.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Geometry and essential chemical and mechanical boundary conditions of the single-layer paint system, as taken from Eumelen et al. (2021).The grey area represents the paint layer and the three red circles indicate the metal soap nuclei present in the initial paint configuration.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8
depicts the deformed paint layer (with the deformation plotted at true scale) for different values of the dimensionless time ,

Fig. 8 .
Fig. 8. Single-layer paint system: Time evolution (using the dimensionless time t) of the paint configuration (with the deformation plotted at true scale), illustrating the growth of the metal soap crystals and the discrete crack patterns forming in the paint material.The contour plot variable  refers to the volume fraction of crystalline metal soap, with the lower bound  = 0 designating the paint material (blue colour) and the upper bound  = 1 indicating the metal soap crystals (red colour).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11 .
Fig.11.Multi-layer paint system: Time evolution (using the dimensionless time t) of the paint configuration (with the deformation plotted at true scale), illustrating the growth of the lead soap crystal and the discrete cracks and delamination processes forming in the paint system.The contour plot variable  refers to the volume fraction of crystalline lead soap.The black lines indicate the intact interface between the ground and the paint layers and the white lines designate two symmetrical, delaminating cracks at the interface.The grey area above the interface indicates layer 2, in which no lead ions are present and therefore does not contribute to the formation of metal soap.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. Multi-layer paint system: Time evolution (using the dimensionless time t) of (a) the normalised crack opening  1 of the mode I surface crack, and (b) the normalised delamination length  d of the two mode II delaminating cracks developing along the layer interface.

Table 2
Mechanical properties used in the single-layer paint model illustrated in Fig.7.

Table 3
Mechanical properties used in the multi-layer paint model illustrated in Fig.10.