Molecular dynamics (MD) simulations are highly attractive for studying the influ- ence of interfacial effects, such as the enrichment of components, on the mass transfer through the interface. In a recent work, we have presented a steady-state MD simu- lation method for investigating this phenomenon and tested it using model mixtures with and without interfacial enrichment. The present study extends this work by in- troducing a non-stationary MD simulation method. A rectangular simulation box that contains a mixture of two components 1 + 2 with a vapor phase in the middle and two liquid phases on both sides is used. Starting from a vapor-liquid equilibrium state, a non-stationary molar flux of component 2 is induced by inserting particles of compo- nent 2 into the center of the vapor phase in a pulse-like manner. During the isothermal 1 relaxation process, particles of component 2 pass through the vapor phase, cross the vapor-liquid interface, and enter the liquid phase. The system thereby relaxes into a new vapor-liquid equilibrium state. During the relaxation process, spatially resolved responses for the component densities, the fluxes, and the pressure are sampled. To reduce the noise and provide measures for the uncertainty of the observables, a set of replicas of simulations is carried out. The new simulation method was applied to study mass transfer in two binary Lennard-Jones mixtures: one that exhibits a strong enrichment of the low-boiling component 2 at the vapor-liquid interface and one that shows no enrichment. Even though both mixtures have similar transport coefficients in the bulk phases, the results for the mass transfer differ significantly, indicating that the interfacial enrichment influences the mass transfer. 2 Introduction Mass transfer through vapor-liquid interfaces is a highly important phenomenon: it is ubiq- uitous in nature and applied in many technical processes, including fluid separation processes such as absorption and distillation. While macroscopic mass transfer theories generally con- sider the interface to be a two-dimensional object, it is well known from molecular studies that an interfacial region exists in which the intensive properties change from the value in one bulk phase to that in the other in a continuous manner, involving steep gradients, as the interfacial region is only a few nanometers wide. Interesting phenomena occur in the interfacial region, such as the enrichment of components.1–14The vapor-liquid interfacial re- gion has been studied extensively in equilibrated systems, but only little is known about the influence of mass transfer in that region and, vice versa, on the influence of that region on mass transfer. That influence is usually neglected in macroscopic mass transfer theories, but recently, doubts on the validity of that assumption have been raised,7,10,15–21which makes such studies all the more important. Molecular dynamics (MD) simulations are ideally suited for investigating mass transfer through interfaces on the molecular level. Given the importance of this phenomenon, it is astonishing that, up to recently, no MD simulation methods for studying mass transfer through vapor-liquid interfaces driven by a gradient in the chemical potential have been described in the literature. We have therefore devised and tested two such methods. The first one is based on stationary simulation and has been published recently.22The second one, a non-stationary method, is described in the present work and tested using the same mixtures as in the tests of the first method. In these simulations, we address the question of whether an interfacial enrichment has an influence on the mass transfer rates. The approach we take for this is straightforward: we compare results for fluxes though interfaces in a model mixture with a high interfacial enrichment with those obtained for a similar mixture that has no enrichment. Heat transfer through vapor-liquid interfaces during evaporation and condensation has 3 been frequently investigated in the literature using molecular simulations23–42and contin- uum models such as density functional theory (DFT).31,36,40,42–46This inherently involves mass transfer, but the focus of the studies is different from ours; furthermore, in most of these studies, only pure components were considered. Stationary mass transfer through pores, membranes, crystals ,47–53and homogeneous fluid phases54–56induced by a gradient of the chemical potential has been the subject of several studies. However, these studies do not involve vapor-liquid interfaces. In these studies, the gradient was usually created using the so-called dual-control volume (DCV) method.49,52,57–59In a recent work by our group,22which was already mentioned above, we have adapted the DCV method to induce a stationary flux of a component through the vapor-liquid interface of a binary mixture. The simulation method for quasi-stationary mass transfer through a vapor-liquid interface driven by a gradient in the chemical potential was adapted in a recent study by Rauscher et al.60with the extension of a temperature gradient to investigate mass and heat transfer at interfaces of similar mixtures of simple fluids. Non-stationarymasstransferclosetoandintheinterfacialregionhasbeeninvestigatedin several studies using different techniques: Baidakov et al.61,62have investigated the dynamic formationofvapor-liquidinterfacesinpuresubstancesandbinarymixturesandtheevolution of interfacial properties in the process. In their studies, first a system containing a coexisting vapor and liquid phase was equilibrated. After the equilibration, all particles in the interface layer were instantaneously removed and the non-stationary relaxation process and formation of a new vapor-liquid interface were observed. Bucior et al.63studied evaporation of a pure liquid in a slit pore into a spontaneously created vacuum by moving the pore walls. Nagl et al.15studied the non-stationary mass transfer through liquid-liquid interfaces of different real fluid mixtures using a combination of experiments and theoretical calculations. Braga et al. investigated self-diffusion at liquid-liquid interfaces64and the free energy barrier at vapor-liquidinterfaces65bymovingaprobeparticlethroughtherespectiveinterface. Garrett et al.66observed the trajectories of single particles moving in the vicinity of a vapor-liquid 4 interface in aqueous systems under equilibrium conditions and categorized four different behaviors of these particle: adsorption, absorption, desorption, and reflection. In the present work, we describe a new simulation method for studying a non-stationary mass transfer through vapor-liquid interfaces of mixtures: First, the system containing only particles of one component is equilibrated. The initial equilibrium state consists of a vapor phase slab surrounded by a liquid bulk phase, which coexists in a rectangular simulation box with periodic boundary conditions. After the first equilibration, particles of a second component are inserted in the middle of the vapor phase during a short period of time, i.e., in a pulse-like manner. This perturbation induces a non-stationary, one-dimensional flux of the inserted particles towards the two vapor-liquid interfaces and through them into the liquid phase on both sides of the vapor phase. The process is basically symmetric with regard to these interfaces. The relaxation process caused by the perturbation is then studied, focusing on the processes in the vicinity of the interfaces. Simulations are run until a second, final equilibrium state is reached. A major problem of non-stationary molecular dynamics simulations is to distinguish an actual physical process from random noise and to provide a measure of uncertainty. Therefore, our simulation study was performed using sets of replicas , i.e., the simulation was repeated several times using different initial velocity distributions. The simulation box and start system size, temperature, and the insertion procedure were the same for all simulations in a set of replicas. The results obtained from a set of replicas were averaged, resulting in a reduction of the noise of the signal; furthermore, a measure for the statistical uncertainties of the observables was obtained. This replica strategy based on varying starting conditions is a well-established method in non-stationary molecular dynamics simulations.65,67–70 Thenewsimulationmethodwasappliedtoinvestigatenon-stationarymasstransferintwo binary Lennard-Jones mixtures with distinctively different phase behavior. Both mixtures (referred to as mixture A and B in the following) have been studied extensively by our group regarding thermodynamic10,16,71–73and transport properties.22,74,75Equilibrium molecular 5 dynamics simulations and density gradient theory consistently yield a strong enrichment of the low-boiling component at the interface for mixture A, while no such enrichment is found in mixture B.10,16,22,71–73The same mixtures have been used in the quasi-stationary approach of our recent study.22 An enrichment of low-boiling components at equilibrium conditions has been reported for many mixtures of real fluids and model fluids;1–14see Ref. 76 for a recent review. It has been discussed in various studies7,10,15–21that the enrichment might have an influence on the diffusive mass transfer through the interface (which would be in line with Fick’s law), but to date there is no proof of this hypothesis. No experimental studies of this are available, nor are they to be expected in the foreseeable future since the enrichment occurs in the interfacial region, which is only a few nanometers wide and fluctuates, which hinders the acquisition of meaningful experimental data. Hence, we have to rely on simulation methods for corresponding studies. Thispaperisstructuredas follows: First, the simulation methodisintroduced. Then, the studied test mixtures are described. Next, the simulation results from one temperature are discussed in detail, comparing the results from both mixtures. Finally, results obtained for several temperatures are compared, which provide insights into the temperature dependency of the relaxation process. Simulation Method Overview The basic concept of the simulation method is illustrated in Figure 1. The simulation box has a rectangular shape and periodic boundary conditions in all directions. At the start, the simulation box contains only particles of component 1. The conditions are chosen so that the box contains two phases: a vapor phase in the middle of the simulation box surrounded by a liquid phase. During the entire simulation run, the number of particles of component 1, 6 the temperature, and the volume of the simulation box are kept constant. After establishing the vapor-liquid equilibrium state of the pure component 1 system, it is perturbed by the insertion of particles of component 2, which are inserted into the control volume (CV) in the middle of the vapor phase (cf. Figure 1). The response of the system to that perturbation is evaluated until the system relaxes into a second vapor-liquid equilibrium state. Except for statistical fluctuations, the process is symmetric with respect to a plane in the middle of the simulation box. Hence, a simulation run consists of the following five phases: 1. Initial equilibration phase ( IniEq) 2. Sampling phase of the first vapor-liquid equilibrium state ( Eq1) 3. Insertion phase ( In) 4. Relaxation phase ( Relax) 5. Sampling phase of the second vapor-liquid equilibrium state ( Eq2) The abbreviations in the parentheses are used in the following for brevity. The first two simulation phases IniEqandEq1are carried out in the NVTensemble. Subsequently, in the Inphase, particles of component 2 are inserted into the center of the vapor bulk domain (cf. Figure 1). The insertion of particles of component 2 inside the control volume is performed using a Monte Carlo algorithm54,57,58that adjusts the number of particles of the component 2 based on a prescribed chemical potential of that component µ2,CV. Since the chemical potential is controlled only in a subvolume of the entire simulation box (cf. Figure 1), the Inphase can be considered to be a quasi- µV Tensemble. With the beginning of the particle insertion, a one-dimensional, non-stationary molar fluxj2of particles of component 2 is established. The molar flux is positive on the right of the control volume and negative on the left of the control volume due to the symmetry of the simulation box setup (cf. Figure 1). During the Relaxphase the system re-equilibrates 7 to the second equilibrium state point, which is sampled over the course of the Eq2phase. TheRelaxandEq2phase are again carried out in the NVTensemble. Definition of Observables and Data Processing The simulation box is discretized in z-direction into 1200 equally sized bins. Since a planar interface is studied, no bin discretization in x- and y-direction is applied, i.e., the simulation scenarioisquasi-one-dimensional. Observables χ(τ, z)aresampledin thesebinsas afunction of the z-direction and time τevery 10000 simulation steps. The studied observables are the pressure tensor p, the density ρi, and the mole composition xiof a component i. The rate of change of the density ∂ρi/∂τ(τ, z), which is needed for the determination of the molar fluxji, is determined for each bin separately by numerical differentiation with respect to the time τ. The position of the liquid film and its interfaces is not fixed in the simulation box and can therefore freely move during the simulation run (it fluctuates and changes slightly due to the absorption of component 2). The positions z50=z(ρ50)of the two interfaces as well as the positions z10=z(ρ10)andz90=z(ρ90)were determined by the z-position, where the local total density corresponds to 10 %, 50 %, and 90 % of the bulk density difference, i.e. ρ10=ρ′′+ 0.1(ρ′−ρ′′) ρ50=ρ′′+ 0.5(ρ′−ρ′′) ρ90=ρ′′+ 0.9(ρ′−ρ′′),(1) where ρ′′indicates the total vapor bulk phase density and ρ′indicates the total liquid bulk phase density. Details on the procedure to determine ρ′′andρ′are given in the Supporting Information. In the following, the data processing procedure is described for the right interface only (cf. Figure 1). The interface on the left is evaluated by applying an equivalent procedure 8 and results for both interfaces are averaged. The molar flux j2of the inserted particles of component 2 is determined by solving the molar balance of each bin per time step. Particles of component 2 may enter or leave a bin over its two boundaries, thus there are two molar fluxj2per bin. The molar flux j2is always the net molar flux here. A positive flux is directed in positive z-direction and, accordingly, a negative flux is directed in negative z-direction. In the following, the index n= 1...Ndenotes bins and the flux at the left bin boundary has the index nwhile the flux at the right bin boundary has the index n+1. The molar balance of a bin nfor particles of component 2 at a time step τcan therefore be written as ∂ρ2 ∂τ(τ, n) =−j2(τ, n+1)−j2(τ, n) ∆zbin, (2) where ∆zbinis the size of the bin and j2is the molar flux of component 2. To obtain N+1 unknown molar fluxes of Nbins, Nmolar balances have to be solved (cf. Eq. (2)) with the help of a boundary condition. In this work, the symmetry boundary condition applies, according to which the net flux in the middle of the liquid slab is zero (cf. Figure 1). The position of the middle of the liquid slab is determined to be the midpoint between the positions z50of the two interfaces in each time step and hence may fluctuate during the simulation. The gradient and molar balances are determined in a fixed reference frame in each time step. Only the position, where the symmetry boundary condition is applied, moves according to the position of the middle of the liquid slab. This is done for each simulation in a set of replicas independently and results are averaged at the position of the middle of the liquid slab. The density, composition, and pressure tensor sampled in a bin nat the same time τare averaged over all simulations in a set of replicas. This procedure significantly reduces noise and enables the estimation of the statistical uncertainty for each sampled observable χ. The start and end time of the Inphase are the same for all simulations. The number of inserted particles in a set of replicas varies due to the probabilistic nature of the Monte 9 Carlo algorithm. The number for µ2,CVwas chosen so that the mean number of inserted particles in a set of replicas was 1,200 particles. The distribution of the number of inserted particles in a set of replicas is quite narrow (standard deviation about 1.5 %). Preliminary tests showed that this yields a satisfactory signal-to-noise ratio for the observables of interest such as density, pressure, and molar flux during the simulation run. Details are given in the Supporting Information. The behavior of the system in the vapor phase and the liquid phase in the vicinity of the vapor-liquid interface is of particular interest. The relaxation process is sampled in two measurement volumes (MV) close to the interface: one in the bulk vapor phase and one in the bulk liquid phase (here labeled with MV vapand MV liq, respectively). The positions of these measurement volumes are defined with respect to the interface position, as illustrated in Figure 2. Figure 3 illustrates the simulation procedure, the responses to the perturbation in vapor and liquid phase determined from a set of replicas as well as the five phases of the non-stationary simulation method. Application of the Simulation Method Binary Lennard-Jones Mixtures The Lennard-Jones (LJ) fluid is computationally inexpensive, but is able to reliably describe properties of simple fluids,77,78which makes it well-suited for testing new methods in molec- ular simulations. In the present work, the so-called truncated and shifted Lennard-Jones potential (LJTS) was used, with a cut off radius of rc= 2.5σ: uLJ(r) = 4 εσ r12 −σ r6 and (3) uLJTS(r) =  uLJ(r)−uLJ(rc) r≤2.5σ 0 r >2.5σ. (4) 10 Here, uLJ(r)isthepotentialenergyoftheLJpotentialbetweentheparticles,79whichdepends ontheirdistance r. Forthecross-interactionsofunlikeparticlesinthemixtures, themodified Lorentz-Berthelot combination rules were applied80,81 σij=σi+σj 2(5) εij=ξij√εiεj (6) in which an additional binary parameter ξijis introduced (cf. Eq. (6)), which has an impor- tant influence on the phase behavior of binary mixtures.71 The simulation method proposed in this work is applied to two different binary LJTS mixtures. For both mixtures, the size parameter σiand mass miof component 1 and 2 are equal and set to unity. Also, the energy parameter of component 1 ε1is set to unity while ε2of the components 2 is different for the two mixtures. The same holds true for ξ12. The parameters of the two studied test mixtures are given as •Mixture A: ε2/ε1= 0.6andξ12= 0.85 •Mixture B: ε2/ε1= 0.9andξ12= 1.00 The phase behavior and interfacial properties of these two LJTS mixtures have been investi- gatedsystematicallyinpreviousworksofourgroup.10,16,71–73Theirthermophysicalproperties were therein found to be well described by the PeTS equation of state.10,82The resulting phase equilibria of mixtures A and B at the temperatures that were studied in the present work are shown in Figure 4. Mixture A exhibits an asymmetric, wide-boiling phase behav- ior, whereas mixture B has essentially an ideal phase behavior in the sense of Raoult’s law. Mixture A also exhibits a vapor-liquid-liquid equilibrium (VLLE) at the lowest studied tem- perature. Examples of real mixtures, which exhibit phase behavior qualitatively resembling that of mixture A, are water + carbon dioxide, which is relevant ,for example, for carbon capture technologies, nitrogen + ammonia, which occurs in the Haber-Bosch process, and 11 water + ammonia, which among others is used as a disinfectant or household cleaner. In the case of mixture B, real mixtures that have comparable phase behavior are benzene + toluene and propane + propene, both of which are separated via distillation in the petrochemical industry. Even though both mixtures exhibit distinctively different phase behaviors, they have similar diffusivity (the product of density and the mutual diffusion coefficient) in the vapor and liquid phases.22In mixture A, the low-boiling component exhibits a large en- richment at the vapor-liquid interface in thermodynamic equilibrium,10,72i.e., the number density ρ2(z)shows a distinct maximum at the interface. In mixture B, no such enrichment is present.10,72The two mixtures also differ in other interfacial properties such as the surface tension, the interfacial thickness, and the relative adsorption of the low-boiling component at the interface; for details see Ref. 72. In the following, all physical properties are reduced using the LJ potential parameters of the high-boiling component 1 ( σ1,ε1,m1) (which are equal for both considered mixtures) and the Boltzmann constant kB.79,83 Simulation Details For the simulations carried out in this work, the simulation box had a rectangular shape with a length lzinz-direction of lz= 150 σand a quadratic cross section in x- and y- directions with a side length of ly=lx= 27.85σ. The temperature was controlled by the velocity scaling thermostat.79The control volume had a length in z-direction of ∆zCV=10σ and the same quadratic cross section in x- and y-directions as the simulation box. The time step was ∆τ=0.001σ(m/ε)1/2. The total run time was ∆τSim=4 500 σ(m/ε)1/2for all simulations. Thedurationsofthefivesimulationphaseswas ∆τIniEq= ∆τEq1=500 σ(m/ε)1/2, ∆τIn= 100 σ(m/ε)1/2,∆τRelax= 2 900 σ(m/ε)1/2, and ∆τEq2= 500 σ(m/ε)1/2. During the In phase, 50 MC insertion/removal trials were performed every 40 MD steps within the control volume to prescribe the chemical potential. The vast majority of accepted MC trials were insertion moves, as expected. 12 Eachsetofreplicascontained100simulations. Theinitialconfigurationofasetofreplicas contained 16000 particles of component 1. Observables χof bins at the same position nand timeτof each simulation in a set of replicas were arithmetically averaged and the uncertainty in each observable χof a set of replicas was computed with the single standard deviation. Observables χsampled in each bin were averaged at each ∆τ= 10σ(m/ε)1/2time interval. Each bin had a length in z-direction of ∆zbin=0.125σand the same quadratic cross section inx- and y-directions as the simulation box. Overall, the simulation box was discretized in thez-direction into 1200 bins. The rectangular measurement volumes MV liqand MV vaphad a length in z-direction of ∆zMV= 1σand had a quadratic cross section with the same side length as the simulation box in the x- and y-direction. The two measurement volumes were positioned at a distance of ∆zoff=7σfrom the characteristic interface positions z10andz90 (cf. Eq. (1) and Figure 2), respectively. The molar fluxes in and out of the measurement volumes were sampled at the boundaries of the measurement volume in positive and negative z-direction. In the further course of this work for each measurement volume, only the molar flux at the boundary closest to the interface is shown and discussed. The simulations in this work were carried out using the open source molecular dynamics code ls1.84 Results and Discussion In this work, for each of the two studied mixtures, the same four temperatures were inves- tigated, which yielded 8 sets of replicas, each containing 100 simulations. First, the results for both mixtures for the temperature T=0.715εk−1 Bare presented and discussed in detail. Next, an overview of the results for other temperatures is discussed, details and the corre- sponding numerical data are presented in the Supporting Information. Table 1 summarizes the main simulation settings and some key simulation results for all studied sets of replicas: the temperature Tand the chemical potential µ2,CVprescribed in the control volume as input parameters as well as the mean number of particles N2inserted in the Inphase, the 13 vapor pressure of Eq1andEq2, and the vapor and liquid phase density of component 2 in theEq2phase. Results for the Temperature 0.715εk−1 B Response in the Liquid and the Vapor Phase Figure 5 presents the results obtained in the replica studies carried out for T= 0.715εk−1 B with the two mixtures A and B. Besides the mean value obtained from the averaging of the replica results, the standard deviation is also depicted. Results from the sample volumes on the vapor side (MV vap) and the liquid side (MV liq) are shown for three observables: the density of component 2, the flux of component 2, and the pressure. Let us first discuss the uncertainty of the results. For both mixtures and both measure- ment volumes, the noise of the results for the density is of the order of 10 % (cf. Figure 5), which is low considering the difficulty of the spatially and temporally resolved measurement of a small quantity. Results of a similar quality were obtained for the pressure in the mea- surement volume in the vapor phase. In contrast, the uncertainty in the results for the pressure in the measurement volume in the liquid phase is high, which is not astonishing regarding the extreme sensitivity of this property.85,86As the flux of component 2 was not measured directly but derived from the results for the density, no standard deviation is re- ported. Fast changes in the pressure and the density of component 2, which occur during and shortly after the Inphase, would not be recognizable as physically meaningful without the averaging procedure of the set of replicas. Intheinitialequilibriumstate Eq1,whichisthesameforbothmixtures,onlycomponent1 ispresent. Component2, whichisdifferentformixtureAandB,isaddedonlyintheinsertion phaseIn. Shortly after the beginning of the insertion, the density of component 2 and the pressureinthemeasurementvolumeMV vapincreasesteeply. Thereisashorttimedelayuntil the first component 2 particles reach the measurement volume MV vap. As a consequence, the molar flux of component 2 j2in the measurement volume MV vapalso rises steeply in 14 both simulations. All observables shown in Figure 5 peak shortly after the insertion phase has ended. The following relaxation to the second equilibrium state Eq2is different for mixture A andmixtureB,whichispartiallyaconsequenceofthedifferencesintheequilibriumstate Eq2 for the two mixtures. The solubility of light-boiling component 2 in the liquid component 1 is significantly lower for mixture A than for mixture B.71,73This is caused by the lower value ofε2(causing a higher volatility of component 2) as well as by the lower value of ξ12(less favorable mixed dispersive interactions) in mixture A. As a consequence, the final values of the density ρ2in the liquid phase are much lower for mixture A than for mixture B (cf. Figure 5 (b)) and vice versa for the vapor phase (cf. Figure 5 (a)). The larger amount of component 2 in the vapor phase also causes the higher pressure observed for mixture A in the final equilibrium state Eq2(cf. Figure 5 (e) and (f)). As expected, the pressure is the same in both measurement volumes. The lower solubility of component 2 in mixture A also explains why the integral under the molar flux of component 2 in the measurement volume in the liquid phase MV liqshown in Figure 5 (d) is much smaller than the corresponding integral for mixture B. The time constant of the relaxation process is of the same order of magnitude for both studied mixtures. The relaxation for mixture B is slightly slower, which is attributed to the fact that more of component 2 has to be transferred into the liquid phase before the new equilibrium state is reached. The time-dependent changes in the pressure are basically a consequence of the mass transfer, as the mechanical equilibrium is established very fast. As explained above, there is a maximum in all observables shown in Figure 5 in the vapor phase measurement volume MV vapthat is reached shortly after the end of the insertion phaseIn(cf. Figure 5 (a),(c),(e)). After passing through that maximum, for mixture B, all observables simply decay to their values in the second equilibrium state Eq2. At about τ=1 200 σ(m/ε)1/2, a slight change in the way the decay takes places is noticeable, which is visible as a small hump in the signal for j2(cf. Figure 5 (c)). 15 The corresponding findings for mixture A are completely unexpected: after having passed the maximum, the density ρ2measured in MV vapfor mixture A first decreases and then rises again to approach its value in the equilibrium state Eq2(cf. Figure 5 (a)). This oscillating behavior translates into a net flux of component 2 that is first directed towards the vapor- liquid interface ( j2>0) but then, for a short period in time, goes in the reverse direction (j2<0) to then turn back to the expected direction again ( j2>0) (Figure 5 (c)). The negative flux of j2can be interpreted as a reflection of particles of component 2 at the vapor-liquid interface, i.e., an important part of the particles of component 2 reaching the interface cannot enter the liquid phase. This is related to several effects: first, the low solubility of component 2 in the liquid phase of mixture A, which obviously hinders the uptake of component 2 by the liquid phase, and, second, the enrichment that builds up in the interfacial region in mixture A, as will be discussed below in more detail. The low affinity of the liquid phase for component 2 in mixture A is a consequence of the choice of the molecular parameters (low values of ε2andξ12for mixture A). In the discussion of the decay of the signals for mixture B, we have mentioned a slight change at about τ=1 200 σ(m/ε)1/2. This could be due to an effect similar to the one observed for mixture A, which, however, is much milder for mixture B, as the solubility of component 2 in the liquid phase is much better for that mixture. The fact that particles can be repelled from a vapor-liquid interface has been reported by Garret et al.,66who considered, however, equilibrium conditions in which they tracked individual particles (e.g., OH and O 3radicals at water-air interfaces87,88). To the best of our knowledge, this effect has never before been observed in mass transfer studies. Due to the symmetry of the overall simulation scenario, the repelled particles are once again repelled at the symmetry plane in the middle of the vapor phase, which could explain the small oscillations that were observed mainly in the signal for the flux j2in the vapor phase measurement volume MV vap(cf. Figure 5 (c)). Figure 6 shows the transient response in the pressure-composition diagram at T= 16 0.715εk−1 Bfor both studied mixtures. The pressure sampled in the liquid and vapor mea- surement volumes is shown as a function of composition x2and time τ. The liquid and vapor phase state points of both Eq1(before the insertion) and Eq2(after the insertion and relax- ation) have the same pressure, as expected in a vapor-liquid equilibrium. The Eq1andEq2 points of the vapor and liquid domain agree well with the dew and bubble lines calculated with the PeTS EOS.10,82The pressure sampled in the liquid phase has high noise compared to data from the vapor phase, which is typical for pressure sampled in liquid phases (cf. Figure 5 (e) and (f)). Response at the Vapor-Liquid Interface Figure 7 shows spatial profiles of the density ρ2and the molar flux j2of component 2 in the vicinity of the interface and the neighboring bulk phases for both mixtures at T=0.715εk−1 B. The spatial profiles are discretized in time with an interval of ∆τ= 100 σ(m/ε)1/2. The measured spatial profiles for each time step were superimposed at the interface position z50(τ)to decrease the influence from fluctuations of that position. The Inphase starts at τ= 1 000 σ(m/ε)1/2. The profiles measured before this time are essentially identical and represent the initial equilibrium state Eq1, in which no component 2 is present. We will discuss the results for mixture B first. For that mixture, a temporary density peak of component 2 builds up in the interfacial region (cf. Figure 7 (b)). Within the interfacial region, which is basically determined here by density profiles of component 1, this peak is shifted towards the liquid side of the interface. At about τ= 1 400 σ(m/ε)1/2, this temporary density peak reaches its maximum height and then decays and finally vanishes in the equilibrium density profile of component 2 for the equilibrium state Eq2, which shows no extrema. The temporary density peak of component 2 at the interface for mixture B can be in- terpreted as a kind of jamming. The insertion of the particles of component 2 creates an important flux j2in the vapor phase that is directed toward the interface. There is only 17 little friction in the vapor phase, so the particles of component 2 reach the interface with a high directed velocity. At the interface, they are slowed down by the particles in the liquid. The transport of component 2 in the liquid phase is diffusion-controlled and slower than the vapor phase transport. This leads to a temporary accumulation of component 2 near the interface, which, however, vanishes as the overall transport goes to zero upon approaching the equilibrium state Eq2. For mixture B, the molar flux j2is positive at all times and at all locations as expected (the flux is always directed from the vapor phase to the liquid phase). However, the flux is far from being uniform in the volume shown in Figure 7, which extends by 10σto both sides of the vapor-liquid interface. The flux j2measured on the vapor side shortly after the insertion is very large but decays quickly ( τ <500σ(m/ε)1/2) by almost an order of magnitude of the process. The further decay to zero is slower. In interpreting the j2-flux curves, it has to be considered that a negative slope in the profile shown in Figure 7 indicates that component 2 accumulates over time in the considered bin, a constant profile indicates a steady state, and a positive slope indicates depletion. In general, as expected, depletion is observed on the vapor side and accumulation on the liquid side (disregarding the details related to the build-up and vanishing of the density peak discussed above). The exception is the profile measured at τ=1 100 σ(m/ε)1/2on the vapor side where accumulation is observed because the insertion phase just stopped. Also for mixture A, a peak builds up in the concentration profile of component 2 near the interface shortly after the insertion phase Inis finished. But this peak in mixture A is of different nature than the one in mixture B – it does not vanish and also persists in the equilibrium state Eq2. The peak observed in mixture A is caused by the well-known enrichment of light-boiling components at vapor-liquid equilibrium interfaces,1–14which is knowntobeimportantformixtureAbutnotformixtureB.10,72InmixtureA,theenrichment peak builds up very quickly (in about 200σ(m/ε)1/2) and then undergoes only minor changes during the rest of the process (it diminishes slightly, but the difference between the maximal 18 height and the end height are only on the order of 10 %). In contrast to the peak observed in mixture B, the enrichment peak in mixture A is almost symmetric with respect to the interface, which confirms previous findings.72 The flux j2measured in mixture A shortly after the end of the insertion phase In(at τ= 1 100 σ(m/ε)1/2) is similar to that measured in mixture B, which is not unexpected, as the disturbance is then still limited largely to the vapor phase. However, for longer times, major differences occur. First, the flux j2is much lower for mixture A than for mixture B, which can be interpreted as a consequence of the lower solubility of component 2 in the liquid in mixture A. Furthermore, it can be seen that at τ= 1 200 σ(m/ε)1/2the flux j2is reversed in the vapor phase, i.e., there is a net flow of component 2 from the surface back to the vapor. This effect has already been discussed above and is interpreted as a reflection of particles of component 2 at the interface. It may be associated with the buildup of the interfacial enrichment. The temporary density peak of mixture B at the interface position might also be caused by a higher solubility due to the increased pressure. The bubble line of mixture B has a flat positive slope (cf. Figure 6 right). Hence, an increase in pressure significantly increases the mole fraction of component 2 in the liquid phase of mixture B in an equilibrium state, i.e., the solubility of component 2 in the liquid phase increases. The pressure in the vapor phase of mixture B is higher than the pressure of the second equilibrium state Eq2up to τ= 1 700 σ(m/ε)1/2(cf. Figure 5 (d)) while the density peak builds up until τ=1 400 σ(m/ε)1/2 (cf. Figure 7 (b)). The temporarily increased pressure raises the solubility of component 2, which in turn increases the molar flux across the interface and causes the temporary density peak. In contrast, the slope of the bubble line of mixture A is positive and very steep (cf. Figure 6 left). An increase in pressure increases the mole fraction of the liquid phase of mixture A only slightly in an equilibrium state. Mixture A shows only a slight temporary increased density at the interface (cf. Figure 7 (a)). 19 Results for Various Temperatures The new simulation method was used to study the mass transfer process in mixture A and mixture B not only at the temperature T= 0.715εk−1 Bbut also at the temperatures T= 0.66εk−1 B,0.77εk−1 B, and 0.825εk−1 B. This is of particular interest since the interfacial properties of the mixtures are known to depend strongly on the temperature,10,71e.g., the enrichment of mixture A is known to decrease with increasing temperature. Figure 8 shows the results for the density of component 2 as a function of time sampled in the measurement volumes in the vapor phase (MV vap, top) and liquid phase (MV liq, bottom) for both mixtures A and B, for all studied temperatures. In all cases, the change in temperature leads to quantitative changes, but the qualitative behavior is the same for all temperatures, which is not unexpected. Overall, the influence of the temperature on the results is larger for mixture B, mainly as the solubility of component 2 in the liquid depends on temperature for that mixture, whereas it is hardly temperature-dependent for mixture A (cf. Figure 8 (c)). It is interesting that contrary to the overall trend, the height of the initial peak in the density ρ2measured in the vapor phase (MV vap) in mixture A depends strongly on temperature: the peak is large for low temperatures and almost vanishes at the highest studied temperature. Together with the fact that the solubility of component 2 in the liquid phase is hardly temperature-dependent, this is a strong argument for the influence of the enrichment on this peak, namely, as the enrichment is known to be high at low temperatures and vice versa.10,71This argument is further supported by the findings for mixture B for the density ρ2in MV vap, where the influence of the temperature on the peak height is quite low, despite the considerable influence of the temperature on the solubility. Figure 9 shows the flux of component 2 sampled at different positions in the simulation box as a function of time in replica simulations at various temperatures for mixture A (left) and mixture B (right). Results for three positions are reported: the measurement volume on the vapor side (MV vap), the interfacial plane ( z=z50), and the measurement volume on 20 the liquid side (MV liq). The results confirm that the variation in the temperature leads to quantitative changes, but the qualitative behavior remains the same. Differences in the flux of component 2 caused by a variation of the temperature are mainly limited to the first peak in the flux observed shortly after the insertion phase. In remarkable agreement, in basically all cases shown in Figure 8 the peak is highest for the lowest temperature and its height decreases steadily with increasing temperature. The only exception from this rule is the density signal from the liquid side MV liqfor mixture A, which exhibits only a very small peak. Additional results for the temperature dependency of the transient response are given in the Supporting Information. Conclusions A non-stationary molecular dynamics simulation method is introduced for studying mass transfer through vapor-liquid interfaces. The mass transfer of a certain component through theinterfaceisdrivenbyaconcentrationgradientoftherespectivecomponentthatisinitially built up and then the system relaxes in a new equilibrium state. During that relaxation process, the mass transfer through the vapor-liquid interface is examined in detail. The statistics of the sampling of the observables are amplified by studying sets of simulation replicas, which differ only in their initial velocity distribution. In the proposed simulation method, particles of component 2 are inserted over a short period of time in the middle of the vapor phase of an equilibrated vapor-liquid system con- tainingonlyparticlesofcomponent1. Aftertheinsertionphase, thesystemisre-equilibrated. During this relaxation phase, the inserted particles move from the vapor phase through the vapor-liquid interface and enter the liquid phase. Over the whole simulation run, profiles of the density, composition, and molar flux of each component as well as the pressure are sampled. 21 The simulation method was tested on binary mixtures of simple fluids. Two binary LJTS mixtures with the same high-boiling component and different low-boiling components were investigated. Both mixtures exhibit distinctly different phase behaviors and interfacial prop- erties, such as the enrichment,10,71while having a comparable bulk diffusivity.22Mixture A shows a high enrichment, while mixture B shows none. The number of inserted particles is chosen to be the same for the two studied mixtures. Even though both mixtures have similar bulk phase diffusivities, the results from both mixtures from the new simulation method are distinctly different, which can be attributed to differences in the interfacial properties. The behavior of the molar flux across the interface deviates significantly between the two studied mixtures: For mixture A, an important part of the particles is repelled at the interface. The repulsion of particles at the interface causes a temporarily negative net flux back into the vapor phase. No temporary negative flux is observed for mixture B. For mixture B, on the other hand, a temporary density peak is observed at the beginning of the relaxation process on the liquid side of the interface. This temporary accumulation of component 2 particles at the interface is caused by jamming due to a lower transport resistance at the interface and in the liquid bulk phase compared to that in vapor bulk phase. The different behavior of mixture A and B is attributed to three facts: (i) the mass transfer process is governed by the bulk properties of the second equilibrium state, which in turn is dependent on the number of inserted particles and of course the mixture; (ii) a significant amount of particles of component 2 is repelled at the interface in the case of mixture A; (iii) the enrichment acts as a mass transfer resistance. However, the significance of the last two phenomena has to be evaluated in further studies. The results from this work strongly support the hypothesis that interfacial properties influence mass transfer through fluid interfaces. The resistance to mass transfer at the interface is a dynamic property and can eventually become high enough to temporarily reduce or completely stop flux through the interface and accordingly particles are repelled at the interface. The non-stationary simulation method proposed in this work is a valuable 22 tool for further investigations, e.g., for studying state points in the vicinity of a critical point, ternary mixtures, real substance mixtures as well as liquid-liquid interfaces. Conflicts of Interest The authors declare no competing financial interest. Acknowledgement The authors gratefully acknowledge funding of the present work by the ERC Advanced Grant ENRICO (grant agreement No. 694807) and by the Federal Ministry of Education and Research (BMBF, Germany) 16ME0613 WindHPC. The present work was conducted under the auspices of the Boltzmann-Zuse Society of Computational Molecular Engineering (BZS) and the simulations were carried out on the Regional University Computing Center Kaiserslautern (RHRK) under the grant TUK-MTD as well as the Leibniz Supercomputing Centre (LRZ) under the grant (AMSEL)2(pn56mo). Supporting Information Available The Supporting Information is available free of charge at https://pubs.acs.org/doi/10. 1021/acs.jpcb.2c08752 . The Supporting Information contains details on the statistics of two sets of replicas at T= 0.715εk−1 Band additional information on the determination of the position of the interface and the interfacial thickness. Moreover, it contains the simulation results that were not shown in the main body of this work, namely: pressure, density, and flux sampled in the measurement volumes as a function of time at the temperatures T= 0.66εk−1 B,0.77εk−1 B, and 0.825εk−1 Bas well as the corresponding pressure-composition diagrams. Also, the response of the system at the vapor-liquid interfaces for simulations at the temperatures T=0.66εk−1 B,0.77εk−1 B, and 0.825εk−1 Bfor both mixtures are given in the 23 Supporting Information. 24 Tables & figures Figure 1: Scheme of the non-stationary simulation scenario. The snapshot is taken from a sim- ulation during the insertion phase. Prior to the insertion, the system contains only particles of the high-boiling component 1 (red). Particles of component 2 (blue) are inserted into the control volume (blue shaded area), then spread through the vapor phase, cross the vapor-liquid interface, and enter into the liquid phase until a new equilibrium state is established. Figure 2: Schematic of the geometric configuration of the measurement volumes (MV) located in the vicinity of the interface, showing the total density at the interface as a function of the z- direction ( ). The vapor and liquid measurement volumes (blue shaded areas) are located at a distance of ∆zoff= 7σto the characteristic interface points z10andz90, respectively. The posi- tions z10,z50, and z90() are defined according to Eq. (1). 25 Figure 3: Example of results from the evaluation of a set of replicas of simulation runs used here for illustrating the simulation phases. Top: number of inserted particles of component 2 N2as a function of the simulation time (perturbation of the system); Bottom: density of component 2 in the measurement volumes MV vapand MV liq(cf. Figure 2) as a function of the simulation time (response of the system). The dashed lines indicate the equilibrium states at the end of the relax- ation ( Eq2). /s48/s46/s48 /s48/s46/s50 /s48/s46/s52 /s48/s46/s54 /s48/s46/s56 /s49/s46/s48/s48/s46/s48/s48/s46/s49/s48/s46/s50/s48/s46/s51/s112 /s32/s47/s32 /s101/s115/s45/s51 /s120 /s50/s32/s47/s32/s109/s111/s108/s32/s109/s111/s108/s45/s49/s109/s105/s120/s116/s117/s114/s101/s32/s65 /s48/s46/s48 /s48/s46/s50 /s48/s46/s52 /s48/s46/s54 /s48/s46/s56 /s49/s46/s48/s48/s46/s48/s48/s48/s46/s48/s49/s48/s46/s48/s50/s48/s46/s48/s51/s48/s46/s48/s52/s48/s46/s48/s53 /s120 /s50/s32/s47/s32/s109/s111/s108/s32/s109/s111/s108/s45/s49/s84 /s32/s61/s32/s48/s46/s56/s50/s53/s32 /s101 /s107/s32/s45/s49 /s66 /s84 /s32/s61/s32/s48/s46/s54/s54/s48/s32 /s101 /s107/s32/s45/s49 /s66/s84 /s32/s61/s32/s48/s46/s55/s55/s48/s32 /s101 /s107/s32/s45/s49 /s66 /s84 /s32/s61/s32/s48/s46/s55/s49/s53/s32 /s101 /s107/s32/s45/s49 /s66/s109/s105/s120/s116/s117/s114/s101/s32/s66 Figure 4: Phase diagrams of mixture A (left) and mixture B (right) for all investigated tempera- tures calculated with the PeTS EOS.10,82 26 Table 1: Selected input parameters and results from the NEMD simulations for mixtures A and B carried out in the present work. Input: the temperature Tin the whole simulation box and the chemical potential µ2,CVin the control volume during insertion phase. Results: the number of in- serted particles N2as well as vapor-liquid equilibrium properties for the initial and the final state: Pressure pand partial density of component 2 ρ2in the vapor phase (′′) and the liquid phase (′). The uncertainty is reported in parentheses and refers to the last digits of the value. It was ob- tained from the standard deviation from the results from the simulations of the set of replicas. Mixture T µ 2,CV N2 p′′ Eq1 p′′ Eq2 ρ′′ 2,Eq2 ρ′ 2,Eq2 / εk−1 B / ε / εσ−3/ εσ−3/ σ−3/ σ−3 A 0.66 -0.55 1189(43) 0.0030(1) 0.0086(3) 0.0088(4) 0.0106(15) 0.715 -0.675 1199(46) 0.0058(2) 0.0122(3) 0.0092(4) 0.0104(11) 0.77 -0.8 1195(44) 0.0105(2) 0.0175(4) 0.0094(4) 0.0105(8) 0.825 -0.9 1204(45) 0.0173(3) 0.0252(4) 0.0097(5) 0.0103(9) B 0.66 -0.8 1194(66) 0.0030(1) 0.0032(1) 0.0007(1) 0.0538(26) 0.715 -0.9 1212(64) 0.0059(2) 0.0063(2) 0.0011(1) 0.0519(26) 0.77 -1.025 1189(65) 0.0105(2) 0.0111(2) 0.0018(2) 0.0485(26) 0.825 -1.15 1184(59) 0.0173(3) 0.0183(4) 0.0026(2) 0.0455(23) 27 Figure 5: Observables sampled in the replica NEMD simulations in the measurement volumes in the vapor phase MV vap(left) and in the liquid phase MV liq(right), respectively, as a function of the simulation time: density of component 2 ρ2, flux of component 2 j2, and pressure p. The tem- perature was T= 0.715εk−1 B. Results for mixtures A are indicated in red; results for mixtures B in black. Solid lines indicate the mean value obtained from the set of replicas, the shaded area in- dicates standard deviation. The standard deviation is only given for the two properties that were sampled directly ( ρ2andp);j2was calculated from ρ2. 28 Figure 6: Pressure-composition diagram at T= 0.715εk−1 Bfor mixtures A (left) and mixtures B (right). Circles indicate state points sampled in the measurement volume MV liq; triangles indicate state points sampled in the measurement volume MV vap. The bottom plots give a detailed view on the state points sampled in the liquid phase. The color scale indicates the simulation time. Each data point represents the mean value obtained from a set of replicas at a given simulation time τ. The white filled symbols indicate pressure and composition in Eq1(green) and in Eq2 (red) phases. The error bars are the standard deviation obtained from a set of replicas. The black line indicates the phase equilibrium computed with the PeTS EOS.10,82 29 /s48/s46/s48/s48/s48/s46/s48/s50/s48/s46/s48/s52/s48/s46/s48/s54/s48/s46/s48/s56 /s49/s48/s48/s48/s50/s48/s48/s48/s51/s48/s48/s48/s52/s48/s48/s48/s116 /s32/s32/s47/s32 /s115 /s40/s109/s47 /s101 /s41/s49/s47/s50/s114 /s105/s32/s47/s32 /s115/s45/s51 /s114 /s49/s43 /s114 /s50 /s114 /s50/s109/s105/s120/s116/s117/s114/s101/s32/s65 /s82/s101/s108/s97/s120 /s69/s113/s50 /s69/s113/s49 /s73/s110/s114 /s49/s43 /s114 /s50/s114 /s50/s109/s105/s120/s116/s117/s114/s101/s32/s66 /s40/s98/s41 /s40/s97/s41 /s45/s49/s48 /s45/s53 /s48 /s53 /s49/s48/s48/s49/s48/s50/s48/s51/s48 /s40/s122 /s32/s45/s32 /s122 /s53/s48/s41/s32/s47/s32 /s115/s106 /s50/s32/s47/s32/s49/s48/s45/s52 /s32/s115/s45/s51 /s40/s109/s47 /s101 /s41/s45/s49/s47/s50 /s49/s44/s49/s48/s48 /s32/s115 /s40/s109/s47 /s101 /s41/s49/s47/s50 /s49/s44/s50/s48/s48 /s32/s115 /s40/s109/s47 /s101 /s41/s49/s47/s50/s40/s100/s41 /s40/s99/s41 /s45/s49/s48 /s45/s53 /s48 /s53 /s49/s48 /s40/s122 /s32/s45/s32 /s122 /s53/s48/s41/s32/s47/s32 /s115/s49/s44/s49/s48/s48 /s32/s115 /s40/s109/s47 /s101 /s41/s49/s47/s50 /s49/s44/s50/s48/s48 /s32/s115 /s40/s109/s47 /s101 /s41/s49/s47/s50Figure 7: Spatial profiles of the observables sampled in the vicinity of the interface and the neighboring bulk phases for mixture A (left) and mixture B (right). Results for the temperature T= 0.715εk−1 Bfrom sets of replicas. Top: density of the low-boiling component 2 ( ) and total density ( ); Bottom: molar flux of component 2 j2. The simulation time τis indicated by the color. The profiles shown were measured at time intervals of ∆τ=100 σ(m/ε)1/2. 30 /s48/s46/s48/s48/s48/s48/s46/s48/s48/s53/s48/s46/s48/s49/s48/s48/s46/s48/s49/s53/s114/s39/s39 /s50 /s32/s32/s47/s32 /s115/s32/s45/s51/s40/s97/s41 /s77/s86 /s108/s105/s113/s77/s86 /s118/s97/s112 /s48/s46/s48/s48/s48/s48/s46/s48/s48/s53/s48/s46/s48/s49/s48/s48/s46/s48/s49/s53 /s40/s98/s41 /s49/s48/s48/s48 /s50/s48/s48/s48 /s51/s48/s48/s48 /s52/s48/s48/s48/s45/s48/s46/s48/s50/s48/s46/s48/s48/s48/s46/s48/s50/s48/s46/s48/s52/s48/s46/s48/s54 /s116 /s32/s32/s47/s32 /s115 /s40/s109/s47 /s101 /s41/s49/s47/s50/s109/s105/s120/s116/s117/s114/s101/s32/s66 /s109/s105/s120/s116/s117/s114/s101/s32/s65 /s84 /s32/s61/s32/s48/s46/s56/s50/s53/s32 /s101 /s107/s32/s45/s49 /s66 /s84 /s32/s61/s32/s48/s46/s55/s55/s48/s32 /s101 /s107/s32/s45/s49 /s66 /s84 /s32/s61/s32/s48/s46/s55/s49/s53/s32 /s101 /s107/s32/s45/s49 /s66 /s84 /s32/s61/s32/s48/s46/s54/s54/s48/s32 /s101 /s107/s32/s45/s49 /s66 /s40/s100/s41 /s82/s101/s108/s97/s120 /s69/s113/s50 /s69/s113/s49 /s73/s110 /s49/s48/s48/s48 /s50/s48/s48/s48 /s51/s48/s48/s48 /s52/s48/s48/s48/s45/s48/s46/s48/s48/s53/s48/s46/s48/s48/s48/s48/s46/s48/s48/s53/s48/s46/s48/s49/s48/s48/s46/s48/s49/s53 /s116 /s32/s32/s47/s32 /s115 /s40/s109/s47 /s101 /s41/s49/s47/s50/s114/s39 /s50 /s32/s32/s47/s32 /s115/s32/s45/s51/s40/s99/s41 /s82/s101/s108/s97/s120 /s69/s113/s50 /s69/s113/s49 /s73/s110Figure 8: Density of component 2 sampled in the measurement volume in the vapor phase MVvap(top) and in the measurement volume in the liquid phase MV liq(bottom) as a function of the simulation time τ. Results for mixture A (left) and mixture B (right) from replica simula- tions at different temperatures are shown, which are indicated by color. 31 /s48/s53/s49/s48/s106 /s50/s32/s47/s32/s49/s48/s45/s52 /s115/s45/s51 /s40/s109/s47 /s101 /s41/s45/s49/s47/s50/s109/s105/s120/s116/s117/s114/s101/s32/s65 /s48/s49/s48/s50/s48/s51/s48/s109/s105/s120/s116/s117/s114/s101/s32/s66 /s45/s49/s48/s48/s49/s48/s50/s48/s51/s48/s52/s48/s106 /s50/s32/s47/s32/s49/s48/s45/s52 /s115/s45/s51 /s40/s109/s47 /s101 /s41/s45/s49/s47/s50 /s45/s49/s48/s48/s49/s48/s50/s48/s51/s48/s52/s48 /s49/s48/s48/s48 /s50/s48/s48/s48 /s51/s48/s48/s48 /s52/s48/s48/s48/s45/s50/s45/s49/s48/s49/s50/s51/s106 /s50/s32/s47/s32/s49/s48/s45/s52 /s115/s45/s51 /s40/s109/s47 /s101 /s41/s45/s49/s47/s50 /s116 /s32/s32/s47/s32 /s115 /s40/s109/s47 /s101 /s41/s49/s47/s50/s84 /s32/s61/s32/s48/s46/s56/s50/s53/s32 /s101 /s107/s32/s45/s49 /s66 /s84 /s32/s61/s32/s48/s46/s54/s54/s48/s32 /s101 /s107/s32/s45/s49 /s66/s84 /s32/s61/s32/s48/s46/s55/s55/s48/s32 /s101 /s107/s32/s45/s49 /s66 /s84 /s32/s61/s32/s48/s46/s55/s49/s53/s32 /s101 /s107/s32/s45/s49 /s66 /s122 /s32/s61 /s32/s122 /s57/s48/s32/s43 /s32/s55 /s115/s122 /s32/s61 /s32/s122 /s49/s48/s32/s45 /s32/s55 /s115 /s122 /s32/s61 /s32/s122 /s53/s48/s40/s97/s41 /s40/s98/s41 /s40/s99/s41 /s40/s100/s41 /s40/s101/s41 /s40/s102/s41 /s82/s101/s108/s97/s120 /s69/s113/s50 /s69/s113/s49 /s73/s110 /s49/s48/s48/s48 /s50/s48/s48/s48 /s51/s48/s48/s48 /s52/s48/s48/s48/s45/s50/s45/s49/s48/s49/s50/s51 /s69/s113/s50 /s116 /s32/s32/s47/s32 /s115 /s40/s109/s47 /s101 /s41/s49/s47/s50/s82/s101/s108/s97/s120 /s69/s113/s49 /s73/s110Figure 9: Molar flux of component 2 sampled at different positions as a function of time in replica simulations at different temperatures for mixture A (left) and mixture B (right). The tem- perature is represented by the color. Top: measurement volume on the vapor side (MV vap); Mid- dle: interface plane z=z50; Bottom: measurement volume on the liquid side (MV liq). 32