Numerical Modelling of Effects of Biphasic Layers of Corrosion Products to the Degradation of Magnesium Metal In Vitro

Magnesium (Mg) is becoming increasingly popular for orthopaedic implant materials. Its mechanical properties are closer to bone than other implant materials, allowing for more natural healing under stresses experienced during recovery. Being biodegradable, it also eliminates the requirement of further surgery to remove the hardware. However, Mg rapidly corrodes in clinically relevant aqueous environments, compromising its use. This problem can be addressed by alloying the Mg, but challenges remain at optimising the properties of the material for clinical use. In this paper, we present a mathematical model to provide a systematic means of quantitatively predicting Mg corrosion in aqueous environments, providing a means of informing standardisation of in vitro investigation of Mg alloy corrosion to determine implant design parameters. The model describes corrosion through reactions with water, to produce magnesium hydroxide Mg(OH)2, and subsequently with carbon dioxide to form magnesium carbonate MgCO3. The corrosion products produce distinct protective layers around the magnesium block that are modelled as porous media. The resulting model of advection–diffusion equations with multiple moving boundaries was solved numerically using asymptotic expansions to deal with singular cases. The model has few free parameters, and it is shown that these can be tuned to predict a full range of corrosion rates, reflecting differences between pure magnesium or magnesium alloys. Data from practicable in vitro experiments can be used to calibrate the model’s free parameters, from which model simulations using in vivo relevant geometries provide a cheap first step in optimising Mg-based implant materials.


Introduction
Magnesium is a biodegradable, lightweight structured metal that is becoming increasingly popular for orthopaedic implants due to its desirable properties. It is an essential element in the human body for providing normal neurological and muscular functions [1,2], and is detected in large amounts in the bone tissue. It is reported in [3][4][5] that magnesium alloys have a Young's modulus and specific density closer to that of the human bone than the frequently used non-biodegradable titanium and stainless steel implants; this eliminates the problem of stress shielding. Furthermore, the nature of a biodegradable implant saves on expenses and risks associated with undergoing a second surgery to remove the hardware. Biodegradable materials that are typically used in the bone implant industry are polymers and ceramics, which are not as sturdy as metal implants and, consequently, their applications are limited [4], thus highlighting the advantages of a biodegradable metal implant, such as magnesium.
While there are numerous benefits of magnesium implants, in its pure form, magnesium corrodes rapidly in an aqueous environment, which is consistently an obstacle for biomaterial scientists [6]. This corrosion is due to magnesium being reactive to water to form magnesium hydroxide, releasing bubbles of hydrogen that can accumulate to form gas pockets near the implant [4]. Furthermore, the rapid corrosion causes the loss of mechanical support before the newly formed bone tissue can bear the load, thereby preventing the bone from healing correctly. These limitations can be mitigated using alloys of magnesium (e.g., using calcium, zinc, rare earths, etc. [7][8][9]), and research into this for the use in orthopaedic implant devices is vastly growing [10][11][12]. Physiological consequences of magnesium implants has also been explored, including the blood and organ compositions of corrosion products [13] and collagen interaction with implants [14]. However, the major challenge that continues to be faced by regulation and industry is how to predict the corrosion rate of magnesium metal-based biomaterials in vitro and correlate to the timing of its absorption in vivo. In this work, the problem is approached using mathematical modelling aimed to provide a systematic means of quantitatively describing the physiochemical interaction during magnesium corrosion processes in vitro, further informing standardisation of in vitro investigation of magnesium alloy corrosion and implant design parameters for optimal bone growth.
The application of mathematical modelling in metal corrosion has been studied widely for some time [15][16][17][18]. Our approach employs a porous media extension of a model that is used to describe atmospheric corrosion of a block of copper [16,19]; their model is built upon the chemical reactions incurred as the copper sample corrodes using an advection-diffusion system. There are numerous studies of similar problems for a range of metals that neglect the advective contribution to the corrosion dynamics, resulting in Stefan-like problems to describe corrosion in, for example, zirconium [20] and steel [21]. Magnesium corrosion has been the subject of a small number of modelling studies [22][23][24]. In [23], a simple two-phase bulk model of magnesium corrosion was proposed and parameter fitted to experimental data; however, the model is not explicit in the products of corrosion. A spatially explicit model of a galvanised magnesium was proposed in [24]; here, the magnesium block was a fixed domain and they showed that the thickness of the galvanised layer affected corrosion. The authors of [22] used a level-set approach to describe moving interfaces separating a pure magnesium block and a partially corroded phase consisting of dissolved Mg ions (magnesium chloride) and a protective layer of magnesium hydroxide.
In aqueous and in physiological relevant environments, carbon dioxide and dissolved bicarbonates can react with magnesium hydroxide to ultimately form magnesium carbonate, which is largely insoluble and forms a robust protective layer around the magnesium block. This latter feature is currently absent in magnesium corrosion models and is believed to play an important part in the performance of an implant in vivo. In this paper, we adopted the approach of [16], but applied it to magnesium and the resulting layers of corrosive products Mg(OH) 2 and MgCO 3 ; this is the first study to consider the latter product in a magnesium corrosion model. A further novelty is to consider these layers as porous media, whereby there is fluid phase flow within the pores of the developing crystal structures, so that the reactants H 2 O and CO 2 can advect, as well as diffuse, through them. Whether or not this porous media assumption leads to substantially different results will be one of the aspects explored in this paper. An aim is to guide relatively simple in vitro experimentation that can inform the model parameters, which can then be used in the modelling of magnesium in more clinically relevant environments, with more appropriate geometries and dimensions, to predict corrosion in vivo.
In the next section, a partial differential equation (PDE) model is developed as an advection-diffusion system, whereby magnesium is assumed to corrode through a series of chemical reactions in vitro. The model is simplified and non-dimensionalised in preparation for the numerical investigations described in Section 3. In the final sections, the main results are discussed and concluded.

Mathematical Model
In vitro, either in water or a clinically relevant media, the magnesium component of a proposed implant will initially corrode through the reaction with water, according to leading to the production of hydrogen and a protective layer of magnesium hydroxide. The hydrogen gas evolves from the solution leaving the hydroxide film on the magnesium surface, which is only stable at a high pH and, in physiological environments, is vulnerable to further corrosion [11]. In water, dissolved CO 2 reacts with Mg(OH) 2 to ultimately form magnesium carbonate in a reaction summarised by [25,26], which is used as the basis for the model formulation. In water, the principle reacting agent with Mg(OH) 2 are hydrogen carbonates ions, HCO -3 , formed from the reaction of dissolved CO 2 with water. The reaction of these ions with Mg(OH) 2 leads to the formation of magnesium hydrogen carbonate, Mg(HCO 3 ) 2 , which then decomposes to form magnesium carbonate, MgCO 3 . The intermediate magnesium hydrogen carbonate is thermodynamically unstable at atmospheric levels of CO 2 [26] and demonstrated experimentally [7,27]; we thus assume the intermediate hydrogen carbonate form is short-lived and will therefore be neglected in the modelling. A more realistic representation of the overall reaction is [7], In more clinically relevant environments, for example using cell culture medium in vitro, the presence of a hydrogen carbonate buffering system using concentrations reflecting that in the blood (27 mmol/L) [28] means that HCO -3 will lead to the corrosion of Mg(OH) 2 via the reaction in Equation (4). The formulation of the model with regards to conversion of the hydroxide to carbonate forms means that both reactions Equations (2) and (4) are described, and the variable C 2 in this model can be viewed either as the concentration of CO 2 or HCO -3 or both as the stoichiometry for water is the same; for simplicity, the discussion on Mg(OH) 2 corrosion in the remainder of this paper will refer to CO 2 and the reaction in Equation (2). We note that the resulting layer of magnesium carbonate is more stable and has been proposed as a layer to delay the corrosion process [25]. We further assume in the model that, throughout the corrosion process of Mg, the environment is stable (e.g., pH is unchanged, as would be expected in a buffered medium in vitro) and that supplies of water and CO 2 are inexhaustible.
The modelling is aimed at describing events from the beginnings of the corrosion process of a magnesium or magnesium alloy pellet, through the formation of hydroxide and carbonate layers, until the magnesium and magnesium hydroxide is exhausted and only the carbonate remains. For simplicity, we will refer to the magnesium (alloy) pellet as "pure magnesium" to distinguish it from the hydroxide and carbonate forms. The magnesium pellet is assumed to be non-porous and corrodes in a way that maintains a smooth surface, i.e., surface pitting and cracking is assumed negligible at leading order. It is therefore envisioned that, in the corroding process, the magnesium is surrounded by a layer of Mg(OH) 2 , which in turn is surrounded by a layer of MgCO 3 (see Figure 1). In order for the magnesium pellet to corrode further, water must be able to diffuse though the carbonate and hydroxide layers to react at the pellet's surface and CO 2 must be able to diffuse through the carbonate layer to reach the hydroxide compound interface. These assumptions lead to a model that describes both the transport and reaction processes of water and carbon dioxide as well as the location of the interfaces between magnesium and its compounds, which are deposited on the surface of the magnesium metal as corrosion products. The modelling will be formulated for a general 1D geometry, namely Cartesian (describing a magnesium slab), cylindrical (a magnesium rod) and spherical geometry (a magnesium ball). The hydroxide and carbonate layers are treated as porous media, thereby the movement speed, v s i , of the "solid" components, i.e., the Mg(OH) 2 and MgCO 3 , is distinct to that of the fluid and dissolved gas components, i.e., H 2 O and CO 2 , namely v f ; this is a novel feature in metal corrosion models. Fortunately, by assuming ideal geometries, a closed system of equations can be derived based on mass conservation alone.  Figure 1. A physiochemical schematic of the magnesium corrosion system used in the model for cylindrical and spherical geometries. The pure magnesium or magnesium alloy exists in the core (Zone 0, 0 ≤ r < α(t)), the magnesium hydroxide forms a middle layer (Zone 1, α(t) < r < β(t)) and the outer layer consists of magnesium carbonate (Zone 2, β(t) < r < S(t)). Figure 1 shows a cross-section of the cylindrical/spheroid scenarios for the model; the Cartesian case simply has three layers bounded by parallel, linear interfaces. Applying the above ideal geometries means that the resulting model will only consider changes in one-spatial dimension. Writing r and t as the spatial coordinate and time variable, respectively, we denote the coordinates of the moving interfaces as follows:

Mathematical Modelling
• r = α(t) is the location of the magnesium to magnesium hydroxide interface. • r = β(t) is the location of the magnesium hydroxide to magnesium carbonate interface. • r = S(t) is the location of outer edge, exposed to concentrations of water and carbon dioxide, representative of the in vitro environment containing HCO -3 /CO 2 buffering system.
The movement of the magnesium compounds, which is when the hydroxide deposits on the magnesium surface and when the carbonate produces on the hydroxides surface, are denoted by velocities v s i , where i = 1, 2 corresponding to the zone; and the flow of water and carbon dioxide in the fluid phases are denoted by velocities v f i .
It is assumed that the solid structure of Zones 0-2 is homogeneous, i.e., they consist of a fixed volume fraction of the magnesium compound and non-traversable space (ε i ) and traversable space (1 − ε i ); traversable space is defined to be continuous channels of space and excludes completely enclosed gaps in the solid structure. In the magnesium layer, it is assumed that it is entirely non-traversable, hence ε 0 = 1. All reactions are assumed to occur only at the interfaces α(t) and β(t), whereby the rate at which these interfaces move depends on the local rate of reaction. Considering the availability of inexhaustible CO 2 in vitro, we note that in the physiochemical representative corrosion system, the magnesium and magnesium hydroxide regions will eventually vanish, i.e., α = 0 and β = 0, respectively; consequently, there are distinct time phases in the corrosion process that need to be separately handled by the model. We define t = T α as the point in time when α(t) = 0 (i.e., α(t) > 0 for t < T α ), and likewise t = T β for when β(t) = 0. Once β = 0, i.e., for t > T β , there are no further developments in the system and all that remains is a block of magnesium carbonate.
With v s 1 (r, t) and v s 2 (r, t) being the radial velocities of the solid phases, the conservation of mass implies where d = 0, 1, 2 represent Cartesian, cylindrical and spherical geometry, respectively. Here, ε 1 and ε 2 are the solid phase volume fractions, respectively; we note that ε 1 and ε 2 are constant in their respective zones so can be divided out; however, the fraction term will be retained in the model derivation for completeness. We will for simplicity assume that the fluid phase consists of all non-solid materials and that it is non-compressible. Conservation of total material volume implies that for i = 1, 2; hence, using Equations (5) and (6), we have We emphasise that the flow here does not encompass that of the fluid/gas trapped in non-traversable pores in the solid structure.
It is assumed that water can be transported via diffusion and advection throughout Zones 1 and 2, whilst carbon dioxide is limited to Zone 2; carbon dioxide is assumed to be immediately exhausted on contact with Mg(OH) 2 on the r = β interface (see Section 2.1.1). Let W 1 be the mass concentration of water in the pores of Mg(OH) 2 structure and W 2 and C 2 be that of water and carbon dioxide, respectively, in the MgCO 3 . The transport equations for the water and carbon dioxide are with the fluxes defined as J X i = −D X ∂X i /∂r + v f i X i ; again with the ε i s being constants, then the fluid fraction can be divided out in each of the equations. Here, D X is the diffusion coefficient with X representing water, W, or carbon dioxide, C.

Boundary and Interface Conditions
It is assumed that the initial state consists only of a magnesium, and impose where the initial thickness or radius S 0 > 0. Water and carbon dioxide is sourced at the outer surface r = S, which moves at speed equal to the local velocity, hence whereṠ = dS/dt. The conditions on the interfaces are more complex and change at critical points of the corrosion process. On r = α(t), water reacts with the surface of the magnesium block. Assuming that the magnesium surface is uniform, then the rate of reaction, R α , will be dependent on the water concentration and flux there. As two water molecules are consumed, we assume by the law of mass action applied to Equation (1) that R α = kW 2 1 , where k is a rate constant. There are two cases that will be considered in this paper: • Case 1 considers the limit k → ∞, where the reaction is so rapid that water is immediately exhausted on r = α(t) interface, hence W 1 = 0 here. This assumption is most consistent with that used for carbon dioxide on r = β(t). The boundary conditions Equations (13) and (15) are relevant.
1, the small distance for the carbon dioxide to diffuse means that Mg(OH) 2 immediately becomes exhausted on production and Mg converts to MgCO 3 , in effect, immediately; consequently, α(t) = β(t) during this transient. In time, the thickness of Zone 2, S − β, becomes sufficiently large for the reaction to exhaust the carbon dioxide on r = β, allowing the Mg(OH) 2 layer to grow. Let t = T α=β be the smallest time at which C 2 (β(t), t) = 0; then, for t < T α=β , the conditions in Equation (16) hold, whilst, for t > T α=β , Equations (13) and (15) are then imposed.
In both cases, the Mg will eventually be exhausted and the conditions in Equations (14) and (15) are then relevant.
Let A be the area of a surface element on a magnesium surface; then, the volume change rate upon this element of the Mg block is Aα, translating to a molar change rate of µ 0α A (see Table 1). Consequently, the water molar flux through r = α(t) as the boundary moves is (1 − ε 1 )(−W 1α + J W 1 )A/M W = 2µ 0α A, since two molecules of water are consumed and noting constant M W is equal to mass/mol of water. For the k → ∞ case, this provides the equation for the moving boundary α(t), whilst, for k < ∞ ,we have in addition the mass flux through Table 1, the quantities ρ i and µ i represent the values from the magnesium compounds deposited as layers of corrosion products, hence volume elements are inclusive of the void fraction. The volume fraction difference through converting Mg to Mg(OH) 2 is ω α − 1, where ω α = µ 0 /µ 1 ; consequently, volume gain rate from the reaction yields v s 1 A = −(ω α − 1)αA, noting that ω α > 1 implies a gain in volume so that v s 1 must have an opposite sign toα. The final condition results from a no slip condition to the fluid phase on α(t), i.e., v f i =α. In summary, the conditions are on These conditions hold for α(t) > 0. On exhaustion of the pure Mg block and α(t) ≡ 0, the boundary conditions are For the cases on k described above, and t > T α=β for k < ∞, then, on r = β, the carbon dioxide is assumed to be completely consumed by its reaction with Mg(OH) 2 , whilst a water molecule is produced; for the latter, we assume the concentration is continuous across the interface, i.e., Letting A be again the area of a surface element on r = β, then the rate of volume loss of Mg(OH) 2 is (v s 1 −β)A and the molar loss rate is therefore

by definition, and thus the volume gain rate
for β(t) > 0. For k < ∞ and t < T α=β , in which Zone 1 is absent, the boundary conditions are Here, the conversion of Mg to MgCO 3 generates a volume fraction difference of ω β ω α − 1 and the stated condition on v s 1 is equivalent to that in Equation (13) and likewise for v f 1 . The flux condition on water results from the net loss of one molecule from the overall reaction and likewise for CO 2 .
The current model assumes that MgCO 3 will have no exit from the system, so the final state corresponds to when β = 0, whereby all of the Mg and Mg(OH) 2 have been exhausted. Table 2 shows the boundary conditions used in Cases 1 and 2. Note, in the case k < ∞ and T α = T α=β , Phases 2.1 and 2.3 are relevant. Conditions Equations (11) and (12) are relevant for both these cases on k.

Exact Solutions
Equations (5)- (8) are straightforward to integrate, though their solution depends on the various scenarios stated above. For Case 1, applying Equations (13) and (15) yields for t > T α using Equations (14) and (15). For Case 2, applying Equations (15) and (16) gives for t ≤ T α=β , and the velocities for T α=β < t < T α and t > T α are then Equations (17) and (18), respectively. Using the formulation in Equation (17) for v s 2 and the boundary conditionṠ = v s 2 (S, t), we obtain this formula is correct for t < T α=β in Case 2 on insertion of α = β and in the final phase for both cases on substitution of α = 0. From this, we can deduce the final size, S ∞ , of the magnesium carbonate block on substitution of α = β = 0 into Equation (20), giving this can be calculated a priori from the total volume fraction change from converting Mg to MgCO 3 being (S ∞ /S 0 ) 1+d = ω α ω β . There are no further exact solutions to be obtained and variables α, β, W 1 , W 2 and C 2 need to be resolved numerically from Equations (9), (10) and (11)- (16). Table 3 shows the data used for each of the parameters other than initial radius S 0 . There are three parameters ε 1 , ε 2 and k for which information appears to be limited; these will be investigated further for their effects on the degradation behaviour of magnesium metal and the dynamics of biphasic corrosion layers.

Non-Dimensionalisation
Of the various candidates for a suitable scaling in time, none stand out as providing any particular advantage here, we choose the approximate timescale for which carbon dioxide diffuses across a reference distance S * 0 (choosing S * 0 = 1 cm). We rescale the water and carbon dioxide variables with ambient mass concentrations W * 0 and C * 0 and hence write and v * = D Cv * /S * 0 , where the quantities with hats are dimensionless. Using the data in Table 3, the scaling implies thatt = 1 represents about 14.5 h. Let noting that ω α and ω β are already dimensionless; then, on dropping the hats for clarity, we obtain the system with velocities for each scenario the same as Equations (17)- (19). Initial conditions are For Case 1, k → ∞, we have the following boundary conditions for t < T α , interface conditions, and analytical solution for S, and for t > T α , with analytical solution For Case 2 and t ≤ T α=β , we have the following boundary conditions: with interface conditions, and analytical solution For t > T α=β , we have, in addition to Equation (32), the boundary conditions whereby, for T α=β < t < T α , we have Equations (28) and (29) and, for t > T α , we impose Equations (30) and (31). Table 4 displays the values for the dimensionless parameters used in the model simulations in the next section.  Table 3 and Equation (22); † being the free parameters.

Numerical Method and Results
The spacial domains of the system of PDE Equations (23)- (25) are mapped to the unit interval ρ ∈ [1, 2] using the rescaling outlined in Appendix A; the difficulties from the singularities resulting from α = β = S at t = 0 at the start of Phases 1.1 and 2.1 and α = β at t = T α=β at the start of Phase 2.2 are discussed in Appendix B.
The system of PDE Equations (A1)-(A3) and the appropriate boundary conditions for each of the phases was solved using the method of lines [33] implemented in MATLAB (R2017a, MathWorks). The domains for Zones 1 and 2 are divided into a uniform mesh, not necessarily using the same number of points, and the spatial derivatives are discretised using central differences; upwind scheme for the advection terms was also implemented but usually ran slower. The stiff ODE solver, ode15s, was used for the time stepping process.
From the data values listed in Table 3, there is current uncertainty on appropriate values for ε 1 , ε 2 and κ (though ε 1 = 0.6 and ε 2 = 0.4 was chosen for most simulations). We therefore investigate in Sections 3.2-3.5, for the three principle geometries, the effect of these parameters on the model solutions, in particular on the degradation times of the original Mg block and the Mg(OH) 2 layer. We investigate in Section 3.3 the effect of the initial size of the block, which is, of course, an experimentally controllable parameter. In Section 3.6, the significance of the porous media assumption in the current model is examined. We note that all of the results shown are the dimensionless form of the variables, whereby one space unit represents 1 cm and one time unit represents about 14.5 h.

Magnesium Degradation
An example simulation using a finite reaction rate κ (Case 2) is shown in Figure 2 using cylindrical geometry, where κ ≈ 0.04, (corresponding to k = 0.07 cm 4 /g · s in Table 3), ε 1 = 0.6 and ε 2 = 0.4. The size of Mg and its compounds over time are displayed in Figure 2 with the dashed lines showing T α=β , T α and T β . We note that the Mg block degrades relatively quickly at t = O(10), whilst the Mg(OH) 2 takes t = O(10 3 ). This is largely due to a relatively low concentration of CO 2 compared to H 2 O in the fluid phase. Figure 3 displays water and CO 2 concentration distribution at the start of Phases 2.1, end of Phase 2.1 (t = T α=β ), Phase 2.2 at the point the full system is solved numerically (see Section B.3), the start of Phase 2.3 (t = T α ) and the end of Phase 2.3 (t = T β ); the times t are detailed in the caption. In a short time, there is only a very narrow MgCO 3 layer present and, as expected, the water and CO 2 are very nearly uniform r ∈ (β, S). Furthermore, CO 2 is not initially exhausted by the conversion reaction of Mg(OH) 2 to MgCO 3 at r = β, but, in time, it descends, reaching zero as Phase 2.2 begins. As time advances, clear gradients in concentrations emerge and, whilst β − α and S − β remain small, the concentrations of water and CO 2 appear linear. The upward kink in the water distribution is due to production at r = β in the conversion reaction Equation (2), even exceeding the exterior concentration as, locally, water replaces CO 2 molecules. During Phase 2.3, when there is no more Mg remaining, the water distribution W 1 tends to a uniform distribution via diffusion and the zero flux condition at r = 0. By the end of Phase 2.3, the profiles of W 2 and C 2 are no longer linear and CO 2 concentration forms a boundary layer in the vicinity of r = β.  Table 4 and S 0 = 1. The dashed lines show t = T α=β (left), t = T α (middle) and T β (right).    Table 4. Figure 4 shows α, β and S over time for a Case 1 (κ → ∞) example, using fixed solid fractions, ε 1 = 0.6 and ε 2 = 0.4 and an initial radius of S 0 = 1 (representing 1 cm). The results are displayed left to right for Cartesian, cylindrical and spherical geometries. The dashed lines separate the two-phases, Phase 1.1 and Phase 1.2. Here, the geometry is such that the size of the Mg block is greater in the Cartesian case, thus there is more of it to convert to Mg(OH) 2 and ultimately degrade into MgCO 3 , as can be seen from the final sizes, where S

[d]
∞ , (d = 0, 1, 2 for the Cartesian, cylindrical and spherical geometry, respectively), S ∞ . Comparing the cylindrical case with that of Figure 2, we observe that, as expected, the magnesium layer disappears much faster in κ → ∞ case (T α ≈ 1) than for κ ≈ 0.04 (T α ≈ 33); but we note that it does not significantly affect the overall degradation time of Mg(OH) 2 . In reality, the limit κ → ∞ is unlikely to be realistic for pure or mostly pure magnesium, and represents a metal of low purity. However, as degradation of Mg(OH) 2 is independent of κ, there is very little difference between the t = T β values.  Figure 4. Plots of the dimensionsionless variables α, β and S against t, from left to right, Cartesian, cylindrical and spherical geometry using ε 1 = 0.6, ε 2 = 0.4, κ → ∞, the parameters in Table 4 and S 0 = 1. The dashed lines show t = T α (left) and T β (right).

Effect of Magnesium Block Size
The scaling presented in Section 2.3 is such that the initial magnesium block size of S 0 = 1, used in the simulations up to now, represents 1 cm. Figure 5 plots the relationship between the initial Mg radius and the key degradation timescales T α and T β for each of the three principle geometries in a finite κ case. As expected, they show that these timescales increase with the initial radius of the Mg block, S 0 . The plots suggests a power law relationship of T β ∝ S 2 0 . This can be justified due to the concentration of CO 2 being relatively small, so that 1 γ 2 , and hence the decay of the Mg(OH) 2 is slow. Here,β 1 and hence v f 2 1, and then Equation (25) would be expected to tend to the quasi-steady profile, with ∂ r r d ∂ r C 2 ∼ 0, hence ∂ r C 2 (β, t) ∼ A/β d for some constant A > 0. Writing r = S 0 r and β = S 0 β, then, from Equation (30), we obtain where constant B = β d+1 0 /A(d + 1) > 0, thus T β ∝ S 2 0 as shown in Figure 5.  Figure 5. Plots of the dimensionsionless T α and T β against the initial magnesium block size, S 0 , for each of the principle geometries. The parameters used are ε 1 = 0.6, ε 2 = 0.4, κ ≈ 0.3 (k = 0.5 cm 4 /g day) and parameters in Table 4.

Effect of Porosity of the Mg(OH) 2 and MgCO 3 Layers
In Figure 6, the effect of ε 1 (left) and ε 2 (right) on the time scales T α and T β is shown for κ = 0.3, 6 and κ → ∞. The left side of Figure 6 shows that T α increases with the solid fraction, rising sharply as ε 1 → 1. This is to be expected because as ε 1 increases there is less space for water to flow through the Mg(OH) 2 layer, hence decreasing the rate at which water reaches the Mg interface. However, the solid fraction does not have an impact on the degradation time for Mg(OH) 2 as the transport of CO 2 in the MgCO 3 layer governs this process. The figure emphasises the importance of the hydroxide layer at slowing the degradation of the metal core by impeding the passage of water. The nonlinear relationship as predicted by the current model is in accordance with experimental data [34], whilst decreasing in porosity having the effect of enhancing longevity is consistent with Sun et al. [35].
The right side of Figure 6 shows the effects of changing ε 2 whilst keeping ε 1 fixed at 0.6. The solid fraction of MgCO 3 does not appear to have an effect on T α , but does affect T β . Here, for t < T α , the thickness of the MgCO 3 layer, S − β, is fairly small and appears not to be sufficient to impede significantly the passage of water across it, thus T α remains approximately constant. Although, T α varies between the κ values, it is small compared to T β , and the conversion of Mg(OH) 2 to MgCO 3 being independent of κ means the plots are superimposed. As ε 2 increases, it is CO 2 that is impeded by the smaller void fraction, leading to the sharp rise in time T β as ε 2 → 1. Using the argument in Section 3.3 in formulating Equation (36) for large γ 2 , we then expect T β ∝ 1/(1 − ε 2 ); this relationship matches the numerics very well, suggesting that T β → ∞ as ε 2 → 1 − .  Figure 6. Plots of the dimensionsionless T α and T β against ε 1 (left) and ε 2 (right) for κ = 0.03, 6 and κ → ∞, with the remaining parameters listed in Table 4 and S 0 = 1. Figure 7 displays contours of T α for a range of values for κ and ε 1 whilst keeping ε 2 constant at 0.4 (top plots), and a range of values for κ and ε 2 whilst keeping ε 1 constant at 0.6 (bottom plot). The variation in the parameter κ reflects the different degradation rates across various magnesium alloys. As can be observed from the top plot of Figure 7, the longevity of Mg increases as κ and the void fraction 1 − ε 1 decreases. In agreement with Figure 6, the lower plot shows that changes in ε 2 do not have a significant impact on the degradation time of Mg, but a smaller κ lengthens the degradation time of Mg. Figure 8 displays the effects on T α as κ changes in all three geometries (ε 1 = 0.6 and ε 2 = 0.4). As expected, the time for the total degradation of pure magnesium increases as the reaction rate diminishes, whilst the curves tend to the Case 1 solutions as κ → ∞. For the reasons outlined in Section 3.2, Mg takes the longest time to degrade in Cartesian, then cylindrical and then spherical geometry, though the differences are less noticeable on a logged axis as κ → 0.   Table 4 and S 0 = 1.  Table 4.

Role of Advection
The porous media assumption of the corrosion by-products is a novel feature of the current work in metal corrosion studies. The formation of Mg(OH) 2 and MgCO 3 crystal structures allows transport of water and carbon dioxide through its pores. The separate treatment of the resulting fluid and solid phase velocities is in contrast to [16], in which they assumed that the transport of the diffusive species is supplemented by that of the solid phase motion; this presumably reflects these molecules being somehow connected to and dragged along by the crystal structure. Figure 9 compares the evolution of α, β and S from three choices of advective flux velocities, V i , of water and CO 2 , namely Case (i).
The current model based on porous media assumption (V i = v f i , solid lines). Case (ii).
Zero advective transport (V i = 0, dotted lines), i.e., v f i set to zero in Equations (23)- (25) and in the boundary conditions. Case (iii). Advective transport equal to the solid phase velocity, as in [16] (V i = v s i , dashed lines), i.e., v f i swapped with v s i in Equations (23)- (25) and in the boundary conditions.
The plots show that there is little difference qualitatively with the results between the cases, and the only visible difference being in the stages up to about t = T α . Mg is predicted to degrade slightly faster using the current model's assumptions (V i = v f i ) than the zero advection case (V i = 0) and, in turn, is predicted to be faster than that using V i = v s i as the advection flux. This is due to the signs of the advective fluxes being v f i < 0 < v s i , as the reactions generates local solid volume increases, since w α , w β > 1. Thus, negative advective flux in case (i) implies that there is a background inward drift of the reactants towards the reaction sites and hence the overall corrosion rates will be predicted to be faster than of case (ii), with zero drift, and case (iii), where the reactants are drawn away by the drift. For t > T α , the differences are maintained, but, from the figure, the choice of advective flux appears to have little effect on β over a longer period of time.  [16], dashed). The left plot shows the full evolution of the interfaces and the right plot shows the results around t = T α . The results are using cylindrical geometry, with ε 1 = 0.6, ε 2 = 0.4 and parameters in Table 4 and S 0 = 1.

Discussion
In this paper, a partial-differential equation system, with moving boundaries and interfaces, was used to describe the degradation of a magnesium block in aqueous media. This is a first step of modelling degradation of Mg or Mg alloy based orthopaedic implants in biologically relevant environments. The model considers the diffusion and advection of reactants through porous media in the crystal structures generated by reactions with magnesium and its products. Novel features in terms of metal corrosion modelling is the consideration of porous media flow in the crystal structures and the explicit consideration of MgCO 3 . The model was analysed numerically, but small time asymptotic solutions were needed to deal with singularities at initial and certain time points. In principle, the modelling approach is generic and can be used or adapted to model the corrosive process of any metals or alloys.
The porous media assumption leads to the explicit consideration of the solid phase flow, through manufacturing of crystals at the reaction interfaces, and a fluid phase flow to avoid a vacuum. In 1D, using the classical Cartesian, cylindrical or spherical geometries, a closed system of equations can be derived using mass conservation alone. The flow velocities are analytically solvable, and the advective flow of the reactants in the fluid phase does not complicate the model significantly, from one that assumes diffusion as the only means of reactant transport. The investigation in Section 3.6 compared results using three different advection assumptions, whereby V i = v f i (current model) and V i = v s i (used in [16]), representing extreme cases of physically relevant advective velocity V i . The results showed that the choice of advection term can notably affect the predicted time of pure magnesium degradation, though, in the long term, there is little difference in the predicted results. The results of Section 3.2 show that geometry has a significant effect on timescales for degradation and this can impact the shapes of materials used in an implant. For more accurate predictions of corrosion of complex shapes, the problem needs to be studied in two or three dimensions, as has been considered using different modelling approaches in, for example, [22,35]. Extending the current modelling approach to consider 2D or 3D, non-simple geometries poses non-trivial modelling challenges. The increase in number of variables will require constitutive assumptions on the mechanical properties of the Mg(OH) 2 and MgCO 3 layers to close the system [36].
The model consists of three parameters k, ε 1 and ε 2 that are not readily available from the literature. The results of Sections 3.4 and 3.5 show that these can be tuned to predict a wide range of results in terms of timescales for the vanishing of pure magnesium, t = T α , and Mg(OH) 2 layer, t = T β . However, there is scope for these parameters to be estimated based on appropriate in vitro data. For example, data from time-course measurements of the proportion of constituents of small, spherical magnesium or magnesium alloy beads, immersed in appropriate media. The use of small beads, say around 0.1-1 mm radius (see Section 3.3), should ensure the experiment to be completed in a practicable and cheap time frame, whilst them being spherical enables direct application of the model to calibrate k, ε 1 and ε 2 with the data. The interpretation of k can be extended to the corrosion rate of different quality of magnesium metal and its alloy. In particular, the presence of impurity and the grain size of micro-structures, depending on the preparation methods used, can have dramatic effects on the corrosion rate together with environmental factors [37][38][39][40]. A further experiment involves the use of computed tomography (CT) images, which enables spatial details of the macroscale crystal structure that can be used to obtain direct measurement of ε 1 and ε 2 . The model, with these tuned or determined parameters, provides a starting point to predict the corrosion properties of much larger magnesium pellets and in any 2D or 3D extension outlined above.
The current model describes the corrosion of a smooth Mg block in an aqueous media. The constituents of the media that an orthopaedic implant will be exposed to is more complex and may have significant effects on the corrosion behaviour. For example, lower pH corrodes the Mg(OH) 2 and MgCO 3 , so that the pure Mg is more exposed to the environment and accelerating its corrosion [11,41]. Furthermore, the tougher outer layer of MgCO 3 will itself be corroded and the resulting magnesium ions will eventually disperse and be excreted by the host; the modelling of this corrosion process leads to a modified boundary condition on r = S. Chloride ions in plasma will also react with Mg(OH) 2 to form MgCl 2 [9]; here, the model can be extended to consider two reactive species for Mg(OH) 2 , and assume, for simplicity, that the outer layer consists of an isotropic mixture of MgCO 3 and MgCl 2 . In practice, the magnesium block will be pitted and have holes that will presumably affect its corrosive properties as well [42]; this is currently being explored by the authors.
There is thus plenty of scope to improve the current model, in order to describe more realistically the corrosion properties of magnesium based orthopaedic implants in vivo. Nevertheless, the current model provides a promising initial step into a theoretical understanding of magnesium corrosion, hopefully providing useful insights to help make informed decisions on the experimental direction and design of magnesium based implant materials.

Conclusions
To conclude, we summarise the key points from the preceding sections.

•
Mathematically modelled Mg (or Mg alloy) corrosion in aqueous environments (e.g., cell culture medium) and the corrosion products Mg(OH) 2 and MgCO 3 , forming up to three discrete regions of for Mg, Mg(OH) 2 and MgCO 3 , the boundaries of which move in time.

•
The corrosion products are treated as porous media, whereby fluid (water, CO 2 ) and solid (Mg and its compounds) phases move separately. The reactants are transported via diffusion and advection.

•
The model is an advection-diffusion system with multiple moving boundaries, marking the coordinate of the interfaces between the solid phase species.

•
Many of the parameters are obtainable from literature, leaving three free parameters: the reaction rate between Mg (or Mg alloy) and water (k), and the solid volume fractions of Mg(OH) 2 (ε 1 ) and MgCO 3 (ε 2 ) layers.

•
There are two key timescales for Mg corrosion process, namely that of complete corrosion of the original Mg block (T α ) and Mg(OH) 2 (T β ), so that at T β all that remains is a MgCO 3 block. Numerical solutions demonstrated that, over a wide range of parameters, -T α T β , the original Mg block is short-lived relative to the complete corrosion process. -k and ε 1 affect T α , whilst having a little effect on T β . The latter is affected most by ε 2 .
To substantially prolong Mg presence, Mg alloys must have the effect of reducing the reaction rate k. -geometry has an impact on the corrosion timescales. -complete corrosion is described well by the law T β ∝ S 2 0 /(1 − ε 2 ), where S 0 is the original radius (or size) of the Mg block.

•
Relatively simple in vitro experimentation and CT scans can be used to inform the free parameters, thereby enabling the model to predict the outcome in situations more challenging to undertake experimentally.