2.3. Contact metasomatic - hydrothermal iron ore deposits

2. China's major iron ore types

2.1. Sedimentary - metamorphic iron ore deposits

includes seven types [1–6].

4 Iron Ores and Iron Oxide Materials

40% (w(TFe)).

Huanian in Yunnan Province.

2.2. Late magmatic iron ore deposits

Probably due to large national land area and multi-period orogeny, lots of iron deposits developed in China, and the proven total reserves of iron ores are relatively abundant, but mainly low-grade ores. The research shows that the types of iron ore deposits in China are diverse and the metallogenic epoch spanned from the Archean to the Cenozoic, which mainly

This type of iron deposits is also called banded iron formation (BIF). The mineralization ages were mainly Archean and Paleoproterozoic. The reserves of proven mineral resources were probably more than 33 billion tons, accounting for 55% of the total reserves of the country, which is a very important iron ore in China. According to the ore and gouge mineral association, as well as geological features, it would be subdivided into metamorphic iron-siliceous iron ore deposits and metamorphic carbonate iron ore deposits. Their grades were about 25–

The BIF were affected by regional metamorphism and associated with volcanic-iron-siliceous sedimentary deposits. They mainly developed in the Precambrian metamorphic rocks, mostly large deposits. The Anshan in Liaoning Province, Jidong in Hebei Province, Wutai-Lvliang in Shanxi Province, and central Inner Mongolia are mainly developed area where the iron ore deposits were thought the most representative. The southern margin of Yangtze plate, Qinling

The metamorphic carbonate iron deposits were carbonate-type sedimentary iron deposits suffered to minor regional metamorphism. They mainly developed in the Proterozoic strata, and the southeast Jilin Province is known for yielding this deposit. That led to it also called "Dalizi type iron ore deposit" in China, although it is also developed in Yimen, Eshan,

This type of iron ore deposit is associated with the mafic-ultramafic magmatic intrusions of the iron and its iron minerals are rich in vanadium and titanium, commonly referred to as vanadium-titanium magnetite deposits. The mineralization ages are mainly Paleozoic and Proterozoic, with proven reserves of 9 billion tons, accounting for 15% of the total reserves of the country. According to the metallogensis, they can be subdivided into late magmatic differentiated deposits and late magmatic intrusive iron deposits. Late magmatic differentiated iron deposits are formed by the residual magma rich in iron, vanadium and titanium formed by late differentiation of magma crystals. These deposits are mostly of large in reserve and are mainly distributed in Panzhihua and Xichang area in Sichuan Province. They are often called "Panzhihua-type" iron deposits in China. Late magmatic intrusive iron deposits are latedifferentiated iron-bearing ore fluids that break into or along the rock mass. The reserve of this

orogenic belt, Qilian orogenic belt, and the East Tianshan region were also distributed.

This type of iron ores are associated with the contact metasomatism resulted from intrusive rocks and carbonate rocks, and are formed by the exchange of iron-bearing gas-water solutions. Such deposits generally have typical skarn mineral assemblages and are also known as silicon carbon type deposit. Skarn iron ore ages were mainly Mesozoic, ore grade is generally rich, the reserves are generally small and medium, although there are some large ones. Identified reserves of this type are 8 billion tons, accounting for 13% of the total reserves. Skarn iron ore is widely distributed in China, with such deposits as Handan, Laiwu, Daye and Linfen in the east. In addition, Cuihongshan in Heilongjiang Province, Huanggang in Inner Mongolia, Lizhu in Zhejiang Province, Dading in Guangdong Province, Mulonggou in Xianxi Province and Nixon in Tibet also yield skarn type iron ore deposits.


**1. Introduction**

An ectopic kidney refers to the situation in which developmental arrest of the renal ascent is encountered, leading to the location of this viscus in the pelvic, iliac or abdominal area [1]. Most often, a pelvic kidney is situated opposite to the sacrum bone and below the aortic bifurcation and presents an incomplete rotation [2]. Such cases also present certain modifications with regard to the length and disposition of the ureter and vascular supply: In up to half of these patients, a certain degree of hydronephrosis can be encountered, while the renal arteries can arise from the distal aorta, aortic bifurcation, common or external iliac arteries or even from the inferior mesenteric artery [1,2]. Therefore, its pelvic location makes it susceptible to be invaded by all the tumoral processes of pelvic origin, while the surgical approach can be significantly influenced. The incidence of this anatomical particularity is estimated to be one in 2100 to 3000 cases [3].


*4.2. Linear Interpolation of Meta-Models*

In Section 3, it was proposed to interpolate linearly meta-model predictions to extend the latter to arbitrary anvil velocities. Next, in Figures 10–12, we show a comparison between the predictions of the QoI ∆*R* obtained from meta-model interpolation and full FE simulations.

**Figure 10.** Comparison of meta-models interpolated results and FEM simulations, considering ∆*R* for the Johnson–Cook model.

Observing these plots and the results collected in Table 4, we conclude that the meta-model interpolation for the Johnson–Cook and Zerilli–Armstrong models provides accurate predictions of ∆*R* for arbitrary impact velocities. In contrast, the interpolations of the Arrhenius-type model are not as accurate, possibly due to the relatively higher non-linearity of its constitutive equation, affecting directly the flow stress computation. Without a direct means of verifying this assertion, we might speculate that these non-linearities trigger complex deformation patterns in the anvil once the material enters the plastic regime. However, given that the maximum relative error is below 7 × 10 −2 in all three cases, we can accept the interpolated predictions for the three constitutive models. This choice will result in huge computational savings for the Bayesian calibration.

**Table 4.** Errors in the meta-model predictions of ∆*R* compared with full FE simulations.

**Figure 11.** Comparison of meta-models interpolated results and FEM simulation results, considering ∆*R* for the Zerilli–Armstrong model.

**Figure 12.** Comparison of meta-models interpolated results and FEM simulation results, considering ∆*R* for the Arrhenius-type model.

We have also validated the linear interpolation strategy for the quantity ∆*L*. The results are very similar to the ones obtained for ∆*R* and the interpolation plots are not presented. Table 5 collects the errors made by the meta-model for ∆*L* as compared with the FE solution, leading us to conclude, as for the previous QoI, that the interpolated predictions are accurate enough.

**Table 5.** Errors in the meta-model predictions of ∆*L* compared with full FE simulations.


# **Scientific Advances in STEM From Professor to Students**

Edited by

Yadir Torres Hernández, Manuel Felix Angel, Ana María Beltrán Custodio and Francisco Garcia Moreno

Printed Edition of the Special Issue Published in *Metals, Polymers, Foods, Sustainability and Sensors*

www.mdpi.com/topics


**1. Introduction**

Bacteria, fungi (yeasts and molds), mycobacteria, prions, protozoa, and viruses are common pathogens infecting humans and animals. They typically exist within the host or in the environment. It has been observed that these microorganisms exhibit a notable difference in the natural survivability in the environment, as well as susceptibility to chemical and physical inactivation. For example, under ambient and dried conditions, human coronaviruses seem to lose their infectivity in a matter of several hours to several days [1], whereas endospores and prions may remain infectious for years to decades or even indefinitely [2, 3].

As more and more data have become available regarding the survivability and susceptibility of pathogens to microbicides, it has been observed that the pathogens seem to demonstrate an order of susceptibility to chemical and physical inactivation. E. H. Spaulding first proposed a classification system for the sterilization and disinfection of medical instruments based on the infection risk in 1939 [4]. On the basis of this classification, the concept of a hierarchy of pathogen susceptibility was proposed, in which microorganisms are placed into several groups and ranked from least susceptible to most susceptible. In this hierarchy concept, bacterial spores were ranked the least susceptible, followed by mycobacteria, non-enveloped viruses, fungi, vegetative bacteria, and enveloped viruses. The susceptibility hierarchy was also believed to be related to the biochemical and biophysical characteristics of a pathogen [5, 6].

This hierarchy concept has been slightly modified and expanded over the years. For example, prions were added and considered less susceptible to inactivation by microbicides than bacterial spores; small non-enveloped viruses were considered less susceptible than large non-enveloped viruses; and the order between mycobacteria and small non-enveloped viruses was sometimes reversed (**Figure 1**) [7–10]. Additionally, it has been suggested that the hierarchy concept may be applied either "vertically" (i.e., ranking of susceptibility *between* classes of pathogens) and/or "horizontally" (i.e., ranking of susceptibility *within* a class of pathogens) [11].

The hierarchy concept has been quite useful for enabling scientists to better understand the innate difference among various types of pathogens. In the case of newly emerged pathogens, especially, the hierarchy concept has helped stakeholders design and implement a disinfection strategy swiftly with a reasonable level of confidence. The concept also helps the contaminant control for food, pharmaceutical, and biopharmaceutical products, as it is impractical to test every possible contaminating pathogen, and a robust infectivity assay system may be lacking for certain pathogens (e.g., hepatitis E virus).

Despite its usefulness, the hierarchy concept should be interpreted with caution, as it may oversimply the differences and trending of pathogen susceptibilities. Further examination and refinement of the concept may be necessary; and several important questions should be answered. For example, how often do exceptions to the hierarchy occur and what are the underlying reasons? Could a trending be specific to a given type of chemistry? Is the hierarchy the same between susceptibility to both chemical and physical inactivation? Why do pathogens in the same

#### **Figure 1.**

*Proposed hierarchy of susceptibility of pathogens to microbicides. Note: slightly different versions of the hierarchy concept have been proposed in the literature. Mycobacteria have been placed above small non-enveloped viruses, and molds have been placed above large non-enveloped viruses in certain versions. In some versions, the small and large non-enveloped viruses are combined; and yeasts and molds may be combined.*

*Variability and Relative Order of Susceptibility of Non-Enveloped Viruses to Chemical… DOI: http://dx.doi.org/10.5772/intechopen.102727*

group, or even the same family or genus, sometimes exhibit striking differences in susceptibility? Is there a way to identify and separate reliable/consistent trending versus blurred/variable trending? A deeper look at the efficacy data for various types of microbicidal actives, especially for non-enveloped viruses, may help stakeholders understand the scope, reliability, and limitation of the hierarchy concept so that it can be best utilized.

This chapter reviews the inactivation efficacy data from the literature against non-enveloped viruses for several commonly used types of chemistries, either in formulated or unformulated form, in an effort to generate a separate relative order of susceptibility among these non-enveloped viruses for each type of chemistry and to differentiate consistent versus variable trending. Physical inactivation approaches are not covered in this chapter, although a significant degree of variation also exists for physical treatments. It is not clear that the physical inactivation approaches, in general, are governed by the same hierarchy to susceptibility as is observed for chemical inactivation approaches [12].


*Article* **Prediction of Dispersion Rate of Airborne Nanoparticles in a Gas-Liquid Dual-Microchannel Separated by a Porous Membrane: A Numerical Study**

**Zohreh Sheidaei 1,2, Pooria Akbarzadeh 1,\*, Carlotta Guiducci 2 and Navid Kashaninejad 3,\***

3 Queensland Micro- and Nanotechnology Centre, Nathan Campus, Griffith University, 170 Kessels Road, Brisbane, QLD 4111, Australia

**\*** Correspondence<EMAIL_ADDRESS>or<EMAIL_ADDRESS>(P.A<EMAIL_ADDRESS>(N.K.); Tel.: +98-9122022776 (P.A.); +61-737355391 (N.K.)

**Abstract:** Recently, there has been increasing attention toward inhaled nanoparticles (NPs) to develop inhalation therapies for diseases associated with the pulmonary system and investigate the toxic effects of hazardous environmental particles on human lung health. Taking advantage of microfluidic technology for cell culture applications, lung-on-a-chip devices with great potential in replicating the lung air–blood barrier (ABB) have opened new research insights in preclinical pathology and therapeutic studies associated with aerosol NPs. However, the air interface in such devices has been largely disregarded, leaving a gap in understanding the NPs' dynamics in lung-on-a-chip devices. Here, we develop a numerical parametric study to provide insights into the dynamic behavior of the airborne NPs in a gas–liquid dual-channel lung-on-a-chip device with a porous membrane separating the channels. We develop a finite element multi-physics model to investigate particle tracing in both air and medium phases to replicate the in vivo conditions. Our model considers the impact of fluid flow and geometrical properties on the distribution, deposition, and translocation of NPs with diameters ranging from 10 nm to 900 nm. Our findings suggest that, compared to the aqueous solution of NPs, the aerosol injection of NPs offers more efficient deposition on the substrate of the air channel and higher translocation to the media channel. Comparative studies against accessible data, as well as an experimental study, verify the accuracy of the present numerical analysis. We propose a strategy to optimize the affecting parameters to control the injection and delivery of aerosol particles into the lung-on-chip device depending on the objectives of biomedical investigations and provide optimized values for some specific cases. Therefore, our study can assist scientists and researchers in complementing their experimental investigation in future preclinical studies on pulmonary pathology associated with inhaled hazardous and toxic environmental particles, as well as therapeutic studies for developing inhalation drug delivery.

**Keywords:** lung-on-a-chip; porous membrane; gas–liquid dual-channel chip; nanoparticle; numerical simulation

## **1. Introduction**

As a main component of the lungs, alveoli, with a sizeable capillary interface, enable mass transfer and gas exchange in the human body through their air–blood barrier (ABB) (see Figure 1a) [1]. They are consistently exposed to various inhaled NPs, which due to their small size, can easily transfer throughout the airways and reach them. Some NPs are deposited on the alveolar epithelial cells, while others tend to cross the high-permeable ABB and transfer toward other organs and tissues through the vascular system [2–4]. Over the last few decades, morbidity statistics have demonstrated a significant increase in severe

**Citation:** Sheidaei, Z.; Akbarzadeh, P.; Guiducci, C.; Kashaninejad, N. Prediction of Dispersion Rate of Airborne Nanoparticles in a Gas-Liquid Dual-Microchannel Separated by a Porous Membrane: A Numerical Study. *Micromachines* **2022**, *13*, 2220. https://doi.org/10.3390/ mi13122220

Academic Editor: Pingan Zhu

Received: 21 September 2022 Accepted: 11 December 2022 Published: 14 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

pulmonary disorders, including asthma, emphysema, chronic obstructive pulmonary disease, and lung cancer, and highlighted inhaled exposure of alveoli to ambient toxic NPs as the main cause [5,6]. Moreover, new evidence unexpectedly reveals a growing list of other extrapulmonary diseases associated with the translocation of inhaled toxin NPs from the alveoli to the other organs such as brain, kidney, and liver [7]. On the other hand, this transfer mechanism has been recently taken into account in inhalation therapies to develop inhaled nano-drug delivery systems as a targeted and noninvasive method for the treatment, with fewer side effects, of pulmonary and even other organs' diseases [5]. Therefore, all of this results in a burgeoning demand for pathology, toxicology, and therapeutical studies linked with the inhaled NPs in the alveolar region in order to assess the adverse and beneficial effects of inhaled NPs on the target organ or tissue. vere pulmonary disorders, including asthma, emphysema, chronic obstructive pulmonary disease, and lung cancer, and highlighted inhaled exposure of alveoli to ambient toxic NPs as the main cause [5,6]. Moreover, new evidence unexpectedly reveals a growing list of other extrapulmonary diseases associated with the translocation of inhaled toxin NPs from the alveoli to the other organs such as brain, kidney, and liver [7]. On the other hand, this transfer mechanism has been recently taken into account in inhalation therapies to develop inhaled nano-drug delivery systems as a targeted and noninvasive method for the treatment, with fewer side effects, of pulmonary and even other organs' diseases [5]. Therefore, all of this results in a burgeoning demand for pathology, toxicology, and therapeutical studies linked with the inhaled NPs in the alveolar region in order to assess the adverse and beneficial effects of inhaled NPs on the target organ or tissue.

ABB and transfer toward other organs and tissues through the vascular system [2–4]. Over the last few decades, morbidity statistics have demonstrated a significant increase in se-

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 2 of 18

**Figure 1.** (**a**) Air–blood barrier (ABB) in the alveolar space of human lungs; (**b**) Microfluidic lungon-a-chip device comprises air/blood channels separated by a thin, porous membrane and two side vacuum channels. **Figure 1.** (**a**) Air–blood barrier (ABB) in the alveolar space of human lungs; (**b**) Microfluidic lungon-a-chip device comprises air/blood channels separated by a thin, porous membrane and two side vacuum channels.

Lung-on-a-chip devices have recently drawn attention in biomedical studies owing to their unique potential in bio-mimicking the physiology of the lung alveolar ABB and providing more accuracy in cellular responses as compared to conventional cell-culture models. First introduced by Huh et al. [8] in 2010, the lung-on-a-chip comprises air/blood (alveolar/vascular) microchannels separated by a thin, porous membrane to enable the coculturing of epithelial and endothelial cells neighboring on the opposite sides of the membrane and two side vacuum chambers to emulate a breathing lung (see Figure 1b). Adventitiously, this device has opened new research doors in toxico-pharmacological fields, in addition to promising an enormous capacity in recapitulating complex biological pulmonary malfunctions such as edema formation and thrombosis [9,10]. Lung-on-a-chip devices could be conveniently utilized to inject toxic and medical aerosol NPs to replicate inhaled delivery into the alveoli for pathological and therapeutic studies, avoiding unnecessary clinical trials and animal tests. While the air interface in the lung-on-a-chip makes it tremendously unique among the other common liquid-based cellular microsystems, most of the relevant experimental studies aim at streaming an aqueous solution of NPs [8,11,12], and exposure of cultivated cells to the aerosols has been disregarded. Controlled injection of NPs into the lung-on-a-chip will return appropriate sedimentation and distribution in the target area, which is critical for attaining the desired physiochemical efficacy. However, a big research gap exists in the comprehensive understanding of the dynamics of NPs in the lung-on-a-chip device. Accordingly, the distribution and deposition of particles on the endothelial cells in the air channel, as well as their translocation to the media channel, are imperative subjects that should be examined carefully. It is evident that the dynamic behavior of particles in microdevices is determined by many crucial factors, such as the geometry of the device and properties of fluid flow and particles. Nevertheless, optimizing all these parameters is challenging for an enhanced delivery process. Complications associated with fabricating the lung-on-a-chip device and especially the Lung-on-a-chip devices have recently drawn attention in biomedical studies owing to their unique potential in bio-mimicking the physiology of the lung alveolar ABB and providing more accuracy in cellular responses as compared to conventional cell-culture models. First introduced by Huh et al. [8] in 2010, the lung-on-a-chip comprises air/blood (alveolar/vascular) microchannels separated by a thin, porous membrane to enable the co-culturing of epithelial and endothelial cells neighboring on the opposite sides of the membrane and two side vacuum chambers to emulate a breathing lung (see Figure 1b). Adventitiously, this device has opened new research doors in toxico-pharmacological fields, in addition to promising an enormous capacity in recapitulating complex biological pulmonary malfunctions such as edema formation and thrombosis [9,10]. Lung-on-a-chip devices could be conveniently utilized to inject toxic and medical aerosol NPs to replicate inhaled delivery into the alveoli for pathological and therapeutic studies, avoiding unnecessary clinical trials and animal tests. While the air interface in the lung-on-a-chip makes it tremendously unique among the other common liquid-based cellular microsystems, mostof the relevant experimental studies aim at streaming an aqueous solution of NPs [8,11,12], and exposure of cultivated cells to the aerosols has been disregarded. Controlled injection of NPs into the lung-on-a-chip will return appropriate sedimentation and distribution in the target area, which is critical for attaining the desired physiochemical efficacy. However, a big research gap exists in the comprehensive understanding of the dynamics of NPs in the lung-on-a-chip device. Accordingly, the distribution and deposition of particles on the endothelial cells in the air channel, as well as their translocation to the media channel, are imperative subjects that should be examined carefully. It is evident that the dynamic behavior of particles in microdevices is determined by many crucial factors, such as the geometry of the device and properties of fluid flow and particles. Nevertheless, optimizing all these parameters is challenging for an enhanced delivery process. Complications associated with fabricating the lung-on-a-chip device and especially the membrane with micro-scale pores, aerosol delivery, and transient monitoring of NPs, which consequently lead to costly and

laborious optimization procedures, are the prime suspects responsible for such deficiencies in the field [13–15].

Many efforts have been directed toward employing numerical approaches to overcome challenges associated with experimental analyses of various cell-based microfluidic devices to improve their functionality [16,17]. A few recent investigations are devoted to studying the dynamic behavior of solid particles in the lung-on-a-chip. According to their research objectives, these studies follow two main perspectives, namely, the Eulerian approach, which focuses on the concentration of particles considering the problem governed by typical convection–diffusion phenomena, and the Lagrangian approach, which deals with individual particles, tracing the trajectory of each particle separately. For instance, Frost et al. [18] numerically investigated the effect of a porous membrane, one modeled by molecular diffusivity, on the molecular convection–diffusion transport in a bilayer microfluidic device. Their results showed that the molecular concentration at the outlet of the bottom channel increased by membrane porosity augmentation, while it was invertedly affected by the upper channel height.

Although molecule-sized particles are widely considered in many biomedical studies because of their simple diffusivity, particular attention has also been recently paid to particles with diameters ranging between 10 nm to 10 µm due to some other advantages, such as high uptake capacity by cells [19,20]. Because the Eulerian method investigates the average behavior of particles with less computational cost, it is often favored over its more precise counterpart, i.e., the Lagrangian approach. However, the outcome of these two perspectives would not converge for investigating large particles' dynamics, especially when the total number of released particles is relatively low [21]. Using the Lagrangian approach, Arefi et al. [22] studied the deposition of NPs in an air microchannel under a pulsatile bidirectional flow to simulate breathing patterns. They found that increased airflow rate and breathing frequency raise the deposition rate of particles on the substrate of the channel. In another similar study using finite element simulation, Moghadas et al. [23] reported that airborne delivery of NPs to the cellular region in a microchannel is affected by airflow velocity as well as particle diameter.

The lung-on-a-chip device, which enables particle tracing in air and media phases, allows for simultaneously investigating the particles' deposition and translocation. Thus, it is an essential prerequisite for precise toxicological and pharmaceutical studies. Despite its great potential in mimicking the in vivo conditions, the dynamics of NPs in a typical lung-on-a-chip device with a porous membrane and air/blood channels have never been investigated, to the best of our knowledge. To bridge this research gap, the current paper aims to conduct a comprehensive numerical parametric study to provide a clear insight into the role of various hydrodynamical and geometrical properties on the dynamics of airborne NPs in the lung-on-a-chip system. The impact of fluid-flow velocity, membrane porosity, and particle diameter as dominant factors is examined on the dispersion of the particles in both air and media channels as well as their deposition and transfer rates.

This paper is organized as follows: the numerical modeling of the lung-on-a-chip, the governing equations, and the boundary conditions are introduced in Section 2; comparative studies against accessible data and an experimental study are conducted in Section 3 to establish the accuracy of the present numerical study; a case study and parameter sensitivity are discussed in Section 4; and finally, Section 5 provides the conclusions and summarizes the results.

#### **2. Governing Equations and Boundary Conditions**

In the current work, the dynamics of airborne NPs in the gas/liquid (air/blood) channels of a lung-on-a-chip device separated by a thin, porous membrane are investigated against variations in fluid flow velocity, membrane porosity, and particle diameter. For this purpose, a 2D numerical model is developed using Laminar Flow and Particle Tracing for Fluid Flow modules in COMSOL Multiphysics 5.6 finite element software (see Figure 2) [24] to trace the motion of NPs. For simplicity, the system of governing equations is assumed to

be one-way coupled, considering only the effect of the fluid regime on particles and ignoring the reactive effect due to the tiny size of particles. Therefore, first, the time-dependent velocity profile of the fluids is calculated in the channels, and then the obtained solutions are employed to solve the particle tracing. Here, it should be noted that small dimensions and low fluid velocity generally result in laminar fluid flow in microchannels and also allow the airflow to be treated as an incompressible fluid. Moreover, with a Knudsen number smaller than the critical value (*Kn* ∼= 0.0004 < 0.01), the continuum hypothesis and no-slip velocity condition are valid for the airflow field in the channel [25]. The Knudsen number for the airflow is calculated using the hard-sphere collision model as follows [26]: = 2ℎ = 2√2ଶℎ, (1) where is the mean free path and ℎ is the channel height. In this model, (= 346 pm) demonstrates the collision diameter of fluid molecules, (≅ 290 J. Kିଵ/Kg) and represent the specific gas constant and air density, and is the Boltzmann constant. In many experimental studies, blood flow is mainly replaced by a medium flow to supply the cells with required nutrition in microfluidic devices, which could be conveniently considered as a Newtonian fluid with similar hydrodynamic properties of water at 37 °C [23]. Consequently, the Newtonian and incompressible air and medium phases are

2) [24] to trace the motion of NPs. For simplicity, the system of governing equations is assumed to be one-way coupled, considering only the effect of the fluid regime on particles and ignoring the reactive effect due to the tiny size of particles. Therefore, first, the time-dependent velocity profile of the fluids is calculated in the channels, and then the obtained solutions are employed to solve the particle tracing. Here, it should be noted that small dimensions and low fluid velocity generally result in laminar fluid flow in microchannels and also allow the airflow to be treated as an incompressible fluid. Moreover, with a Knudsen number smaller than the critical value (*Kn* ≅ 0.0004 < 0.01), the continuum hypothesis and no-slip velocity condition are valid for the airflow field in the channel [25]. The Knudsen number for the airflow is calculated using the hard-sphere collision

$$Kn = \frac{\lambda}{2h} = \frac{k\_B}{2\sqrt{2}\pi\sigma^2\rho\_a Rh} \tag{1}$$

where *λ* is the mean free path and *h* is the channel height. In this model, *σ* (= 346 pm) demonstrates the collision diameter of fluid molecules, *R* ∼= 290 J.K−1/Kg and *ρa* represent the specific gas constant and air density, and *kB* is the Boltzmann constant. ቀ డ௨ డ௧ + ∙ ∇ ቁ = −∇ + ∇ଶ + g, (3) where and represent the fluid velocity field and pressure, and is the time, while the subscript ( = a, m) stands for air and media, respectively.

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 4 of 18

model as follows [26]:

mation of drag force (ୢ

**Figure 2.** 2D model of airborne delivery of NPs into a gas-liquid dual-channel microdevice with a thin, porous membrane. **Figure 2.** 2D model of airborne delivery of NPs into a gas-liquid dual-channel microdevice with a thin, porous membrane.

The Lagrangian approach is used to model the dynamics of NPs under the action of gravity and hydrodynamic forces arising from the motion of the fluid: ୮ ௗ௨౦ ௗ௧ = ୢ + ୠ + ൬1 − ఘ౦ ൰ g, (4) where ୮ and ୮ are the particle mass and density, and ୮ indicates the particle veloc-In many experimental studies, blood flow is mainly replaced by a medium flow to supply the cells with required nutrition in microfluidic devices, which could be conveniently considered as a Newtonian fluid with similar hydrodynamic properties of water at 37 ◦C [23]. Consequently, the Newtonian and incompressible air and medium phases are governed by continuity and momentum equations as follows [27]:

$$
\nabla \cdot \boldsymbol{u}\_l = \mathbf{0},\tag{2}
$$

$$
\rho\_i \left( \frac{\partial u\_i}{\partial t} + u\_i \cdot \nabla \, u\_i \right) = -\nabla p\_i + \mu\_i \nabla^2 u\_i + \rho\_i \mathbf{g}\_{\prime} \tag{3}
$$

where *u* and *p* represent the fluid velocity field and pressure, and *t* is the time, while the subscript (*i* = a, m) stands for air and media, respectively.

The Lagrangian approach is used to model the dynamics of NPs under the action of gravity and hydrodynamic forces arising from the motion of the fluid:

$$m\_{\rm P} \frac{du^i\_{\rm P}}{dt} = F^i\_{\rm d} + F^i\_{\rm b} + m\_p \left(1 - \frac{\rho\_i}{\rho\_{\rm P}}\right) \mathbf{g}\_{\rm \prime} \tag{4}$$

where *m*p and *ρ*p are the particle mass and density, and *u i* p indicates the particle velocity in air and media channels, respectively. The right-hand side of Equation (4) is the summation of drag force (*F i* d ), the Brownian force (*F i* b ), and the particle's buoyant weight, while the gravity acceleration is considered to be perpendicular to the flow direction. The drag force in Equation (4) is determined in accordance with Stokes law as follows [28]:

$$F\_{\mathbf{d}}^{i} = \Im \pi \mu\_{i} \left(\mu\_{i} - \mu\_{\mathbf{p}}^{i}\right) d\_{\mathbf{p}\prime} \tag{5}$$

where *d*p is the particle diameter. Also, the Brownian force exerted on particles at each time step taken by the numerical solver could be calculated as [29]:

$$F\_{\mathbf{b}}^{\mathrm{i}} = \zeta\_{\mathrm{i}} \sqrt{6\pi k\_{\mathrm{B}} \mu\_{\mathrm{i}} T\_{\mathrm{i}} d\_{\mathrm{P}} / \Delta t},\tag{6}$$

where *Ti* is the fluid temperature, which is considered to be at a constant level of 37 ◦C, and ∆*t* represents the time step. In this calculation, *ζi* , which indicates the direction of Brownian force, refers to a vector whose components are randomly selected with a Gaussian distribution [29]. Inertial lift force contribution is neglected in Equation (4) since the particle Reynolds number is significantly smaller than 1 (*Re*P = *ρiuid* 2 p/*µi l*), resulting in viscous forces dominancy [30].

Fluid flow boundary conditions in the air and media channels are considered as constant inlet velocity, zero outlet pressure, and no-slip on all solid-fluid interfaces. Furthermore, the top surface of the membrane is employed as the sticky boundary condition for the injected airborne NPs allowing evaluation of particle deposition in the air channel. On the other hand, the pass-through boundary condition is applied to the air-media interfaces to allow NPs to enter the media channel unimpededly. It is worthwhile to note that the numerical simulation is terminated once all the transferred particles flow out of the media channel.

#### **3. Model Verification**

In this section, comparative studies against accessible data as well as an experimental study are conducted to establish the accuracy of the present numerical study. Accordingly, four limiting-case verifications are considered by setting the membrane's porosity equal to zero.

#### *3.1. Verification of Numerical Results with Analytical Data*

In the absence of solid particles, the first verification study compares the dimensionless axial fluid velocity in the middle of the air channel with analytical data for a 2D steady laminar flow between two fixed parallel plates in Ref. [31]. For this purpose, fluid velocity and vertical distance from the substrate (*y*) are normalized with respect to the inlet velocity of 0.3 mm/s and the channel height of 100 µm. As depicted in Figure 3a, a fair agreement is observed between the results of the current study and those of Ref. [31].

#### *3.2. Verification of Numerical Results with Other Numerical Studies*

#### 3.2.1. Verification with a 2D Model

Next, the normalized mean concentration of particles with a diameter of 1 µm along the air channel is compared with that obtained by Saidi et al. [21], who developed a customized numerical code to track 2D dynamics of particles under drag and Brownian forces between two parallel plates. Accordingly, a fully developed inlet flow regime with a mean velocity (*u*) of 0.05 m/s and channel height (*h*) of 1 mm is considered. The concentration of particles is calculated by partitioning the computational space into equal bandwidths along the channel length and then dividing the number of particles in each band to the volume of each bandwidth. Finally, distance from the entrance (*x*) and particle concentration (*C*) are nondimensionalized, respectively, with respect to a specific length of 3*hPe*/16, and initial concentration (*C*0 = 105/m2 ) [21]. *Pe* = 2*uh*/*D* stands for the Peclet number, where *D* = 2.96 × 10−11 m2/s is the diffusion coefficient. The outcome, as shown in Figure 3b, demonstrates a good agreement with Ref. [21].

**Figure 3.** Comparison of the present numerical results with (**a**) analytical data in Ref. [31] for a 2D airflow between two parallel plates. Axial fluid velocity and vertical distance from the substrate are normalized with respect to the inlet velocity of 0.3 mm/s and the channel height of 100 µm. Subsequently, (**b**) numerical data in Ref. [21] is shown for the mean concentration of particles with a diameter of 1 µm along two parallel plates. Distance from the entrance and particle concentration are nondimensionalized, respectively, with respect to a specific length of 3ℎ/16, and initial concentration ( = 10ହ/mଶ ). = 2തℎ/ is the Peclet number, and where (= 2.96 × 10ିଵଵmଶ/s) is the diffusion coefficient. Finally, the (**c**) the numerical result of Ref. [22] for deposition rate of NPs with various diameters on the substrate of an air channel is shown. **Figure 3.** Comparison of the present numerical results with (**a**) analytical data in Ref. [31] for a 2D airflow between two parallel plates. Axial fluid velocity and vertical distance from the substrate are normalized with respect to the inlet velocity of 0.3 mm/s and the channel height of 100 µm. Subsequently, (**b**) numerical data in Ref. [21] is shown for the mean concentration of particles with a diameter of 1 µm along two parallel plates. Distance from the entrance and particle concentration are nondimensionalized, respectively, with respect to a specific length of 3*hPe*/16, and initial concentration (*C*0 = 105/m2 ). *Pe* = 2*uh*/*D* is the Peclet number, and where *D* = 2.96 × 10−11 m2/s is the diffusion coefficient. Finally, the (**c**) the numerical result of Ref. [22] for deposition rate of NPs with various diameters on the substrate of an air channel is shown.

#### *3.2. Verification of Numerical Results with Other Numerical Studies*  3.2.2. Verification with a 3D Model

3.2.1. Verification with a 2D Model Next, the normalized mean concentration of particles with a diameter of 1 µm along the air channel is compared with that obtained by Saidi et al. [21], who developed a customized numerical code to track 2D dynamics of particles under drag and Brownian forces between two parallel plates. Accordingly, a fully developed inlet flow regime with a mean velocity (ത) of 0.05 m/s and channel height (ℎ) of 1 mm is considered. The concentration of particles is calculated by partitioning the computational space into equal bandwidths along the channel length and then dividing the number of particles in each band to the volume of each bandwidth. Finally, distance from the entrance () and particle concentration () are nondimensionalized, respectively, with respect to a specific length of 3ℎ/16, and initial concentration ( = 10ହ/mଶ) [21]. = 2തℎ/ stands for the Peclet Subsequently, the correctness of the developed finite element model is examined to evaluate the deposition of airborne NPs. With the definition of "the ratio between the number of deposited particles on the substrate and the number of total released particles" for the deposition rate (*ε*), the third verification study compares the deposition rate obtained through the current 2D model with that reported for a 3D model developed by Arefi et al. [22]. All adjustments for boundary conditions are taken into account, according to Ref. [22]. Here, the effect of channel depth is also surveyed, enabling the "Shallow Channel" interface in the 2D model. Figure 3c compares the results obtained for a 2D channel with finite and infinite depth and the 3D model of Ref. [22], where the independence of deposition rate and the channel depth as well as an excellent agreement between these three is perceivable.

#### number, where (= 2.96 × 10ିଵଵmଶ/s) is the diffusion coefficient. The outcome, as *3.3. Qualitative Verification of the Numerical Results with Experimental Data*

shown in Figure 3b, demonstrates a good agreement with Ref. [21]. 3.2.2. Verification with a 3D Model Subsequently, the correctness of the developed finite element model is examined to evaluate the deposition of airborne NPs. With the definition of "the ratio between the number of deposited particles on the substrate and the number of total released particles" Besides verification against the analytical and numerical data in the literature, an experimental study was accomplished to support the results of the present numerical model. The experiment was carried out with the injection of particles into a microchannel under a water flow. Therefore, a numerical simulation was developed in accordance with the experimental parametric values to compare the obtained results from the two studies. Accordingly, a polydimethylsiloxane (PDMS) microchannel with rectangular cross-section

was fabricated using a silicone mold (see Figure 4a). The mold was made with maskless ultraviolet (UV) lithography method, as shown in Figure 4b. A silicon wafer coated with a positive photoresist was exposed to direct UV laser to transfer the channel design on it. Subsequently, the wafer was dry etched to engrave the microchannel mold. Then, a 10:1 mixture of PDMS was cast on the mold and the microchannel was sealed with a glass slide to finalize the fabrication procedure (see Figure 4b). A syringe pump was used to inject the particles into the channel under different fluid velocities (see Figure 5 for the experimental setup). In order to have a similar inlet boundary condition to that of the numerical model, the particles were first accumulated at the entrance of the inlet and injected altogether into the channel. Fluid injection was stopped once the particles reached the channel outlet or deposited entirely. was fabricated using a silicone mold (see Figure 4a). The mold was made with maskless ultraviolet (UV) lithography method, as shown in Figure 4b. A silicon wafer coated with a positive photoresist was exposed to direct UV laser to transfer the channel design on it. Subsequently, the wafer was dry etched to engrave the microchannel mold. Then, a 10:1 mixture of PDMS was cast on the mold and the microchannel was sealed with a glass slide to finalize the fabrication procedure (see Figure 4b). A syringe pump was used to inject the particles into the channel under different fluid velocities (see Figure 5 for the experimental setup). In order to have a similar inlet boundary condition to that of the numerical model, the particles were first accumulated at the entrance of the inlet and injected altogether into the channel. Fluid injection was stopped once the particles reached the channel outlet or deposited entirely.

Besides verification against the analytical and numerical data in the literature, an experimental study was accomplished to support the results of the present numerical model. The experiment was carried out with the injection of particles into a microchannel under a water flow. Therefore, a numerical simulation was developed in accordance with the experimental parametric values to compare the obtained results from the two studies. Accordingly, a polydimethylsiloxane (PDMS) microchannel with rectangular cross-section

for the deposition rate (), the third verification study compares the deposition rate obtained through the current 2D model with that reported for a 3D model developed by Arefi et al. [22]. All adjustments for boundary conditions are taken into account, according to Ref. [22]. Here, the effect of channel depth is also surveyed, enabling the "Shallow Channel" interface in the 2D model. Figure 3c compares the results obtained for a 2D channel with finite and infinite depth and the 3D model of Ref. [22], where the independence of deposition rate and the channel depth as well as an excellent agreement between

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 7 of 18

*3.3. Qualitative Verification of the Numerical Results with Experimental Data*

these three is perceivable.

**Figure 4.** (**a**) Microfluidic device, which is bonded to a glass slide and contains a channel with one inlet and one outlet; (**b**) Schematic diagram of the fabrication process of the microfluidic device. **Figure 4.** (**a**) Microfluidic device, which is bonded to a glass slide and contains a channel with one inlet and one outlet; (**b**) Schematic diagram of the fabrication process of the microfluidic device. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 8 of 18

**Figure 5.** Schematic diagram of the experimental setup for injecting the particles into a microfluidic device using a syringe pump and monitoring their dynamics throughout the device with a microscope. **Figure 5.** Schematic diagram of the experimental setup for injecting the particles into a microfluidic device using a syringe pump and monitoring their dynamics throughout the device with a microscope.

Magnetic particles of diameter 1 µm and 10 µm were used in the experiments to investigate their sedimentation on the channel substrate under three different water flow velocities. In addition to more deposition efficiency due to their higher inertia and dominant buoyant weight, better control on injection and easier tracing are other main reasons for selecting these particles. The experimental parameter values, including the physical properties of particles, fluid velocities, and the dimension of the microchannel, are listed Magnetic particles of diameter 1 µm and 10 µm were used in the experiments to investigate their sedimentation on the channel substrate under three different water flow velocities. In addition to more deposition efficiency due to their higher inertia and dominant buoyant weight, better control on injection and easier tracing are other main reasons for selecting these particles. The experimental parameter values, including the physical properties of particles, fluid velocities, and the dimension of the microchannel, are listed in Table 1.

Figure 6 qualitatively compares the trend of particle deposition obtained from the numerical and experimental studies under four different cases. The experimental results are actual images showing the particle deposition rate on the bottom layer of the microchannels. These images were obtained from two different parts of the channel substrate,

**Fluid Velocity (**,**) [mm/s]**

**Flow Rate [µL/min]**

**Dimension of Channel (**××**) [µm3]**

1 2200 2000 × 100 × 15,000 0.03 and 0.1 0.36 and 1.2 10 1470 2000 × 300 × 15,000 0.1 and 0.3 3.6 and 10.8

one from the first half of the channel and the other from its second half.

in Table 1.

**Density of Particles () [kg/m3]**

**Diameter of Particles () [µm]**

**Table 1.** The values of the parameter used for the experimental test.

Figure 6 qualitatively compares the trend of particle deposition obtained from the numerical and experimental studies under four different cases. The experimental results are actual images showing the particle deposition rate on the bottom layer of the microchannels. These images were obtained from two different parts of the channel substrate, one from the first half of the channel and the other from its second half. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 18

**Figure 6.** Qualitative **c**omparison of distribution and deposition of the particles on the microchannel substrate obtained from numerical and experimental studies. (**a**) ୮ = 1 µm, , = 0.03 mm/s; (**b**) ୮ = 1 µm, , = 0.1 mm/s; (**c**) ୮ = 10 µm, , = 0.1 mm/s; and (**d**) ୮ = 10 µm, , = 0.3 mm/s. **Figure 6.** Qualitative comparison of distribution and deposition of the particles on the microchannel substrate obtained from numerical and experimental studies. (**a**) *d*p = 1 µm, *um*,0 = 0.03 mm/s; (**b**) *d*p = 1 µm, *um*,0 = 0.1 mm/s; (**c**) *d*p = 10 µm, *um*,0 = 0.1 mm/s; and (**d**) *d*p = 10 µm, *um*,0 = 0.3 mm/s.

First, the particles with a diameter of 1 µm were injected into the channel with a height of 100 µm under a fluid velocity of 0.3 mm/s. Due to the dominancy of drag force over the buoyant weight, all the particles were washed out from the channel without any deposition on the substrate (data not shown). Then, the dynamics of particles were examined with a lower fluid velocity of 0.1 mm/s and 0.03 mm/s. Under the fluid velocity of 0.1 mm/s, the particles again streamed quickly toward the microchannel outlet, resulting in a low deposition rate along the channel length (see Figure 6a). The fluid velocity, therefore, was reduced to 0.03 mm/s so that the particles could be uniformly distributed all over the channel substrate and demonstrate a higher deposition rate, Figure 6b. First, the particles with a diameter of 1 µm were injected into the channel with a height of 100 µm under a fluid velocity of 0.3 mm/s. Due to the dominancy of drag force over the buoyant weight, all the particles were washed out from the channel without any deposition on the substrate (data not shown). Then, the dynamics of particles were examined with a lower fluid velocity of 0.1 mm/s and 0.03 mm/s. Under the fluid velocity of 0.1 mm/s, the particles again streamed quickly toward the microchannel outlet, resulting in a low deposition rate along the channel length (see Figure 6a). The fluid velocity, therefore, was reduced to 0.03 mm/s so that the particles could be uniformly distributed all over the channel substrate and demonstrate a higher deposition rate, Figure 6b.

There is a mild deposition gradient along the channel width due to the parabolic profile of the fluid velocity in the channel. Moreover, on the substrate next to the channel side walls, the deposition rate is very low as a consequence of the no-slip boundary condition. In the second part of the experiment, particles with a diameter of 10 µm were examined under fluid velocities of 0.1 mm/s and 0.3 mm/s (see Figure 6c,d). The height of the microchannel was increased to 300 µm to allow the particles easily to move inside the channel. The distribution of particles on the substrate of the channel reveals that particles with a diameter of 10 µm under the fluid velocity of 0.1 mm/s provide a concentrated deposition at the channel's entrance. However, an increase in the fluid velocity to 0.3 mm/s drives the particles quickly forward, resulting in a uniform distribution on the channel substrate. The experimental results for deposition and distribution of the particles on the substrate of the microchannel show a good agreement with the numerical results. There is a mild deposition gradient along the channel width due to the parabolic profile of the fluid velocity in the channel. Moreover, on the substrate next to the channel side walls, the deposition rate is very low as a consequence of the no-slip boundary condition. In the second part of the experiment, particles with a diameter of 10 µm were examined under fluid velocities of 0.1 mm/s and 0.3 mm/s (see Figure 6c,d). The height of the microchannel was increased to 300 µm to allow the particles easily to move inside the channel. The distribution of particles on the substrate of the channel reveals that particles with a diameter of 10 µm under the fluid velocity of 0.1 mm/s provide a concentrated deposition at the channel's entrance. However, an increase in the fluid velocity to 0.3 mm/s drives the particles quickly forward, resulting in a uniform distribution on the channel substrate. The experimental results for deposition and distribution of the particles on the substrate of the microchannel show a good agreement with the numerical results.

The current numerical parametric study investigates the dynamics of airborne NPs

Equations (2)–(4). Accordingly, the impact of fluid flow velocity, membrane porosity, and particle diameter, as dominant factors, are analyzed on the dispersion of the particles in both air and media channels as well as their deposition and translocation. Numerical results are confined only to some specific cases due to computational limitations, as well as for the sake of brevity. The diameter of solid particles is assumed to be ranging between 10 to 900 nm, according to their high deposition rate reported in the alveolar region [4]. The effect of membrane porosity is examined for pore diameters () of 3 and 10 µm, and

**4. Numerical Results**

#### **4. Numerical Results**

The current numerical parametric study investigates the dynamics of airborne NPs in a gas–liquid dual-channel lung-on-a-chip device with a thin, porous membrane using Equations (2)–(4). Accordingly, the impact of fluid flow velocity, membrane porosity, and particle diameter, as dominant factors, are analyzed on the dispersion of the particles in both air and media channels as well as their deposition and translocation. Numerical results are confined only to some specific cases due to computational limitations, as well as for the sake of brevity. The diameter of solid particles is assumed to be ranging between 10 to 900 nm, according to their high deposition rate reported in the alveolar region [4]. The effect of membrane porosity is examined for pore diameters (*d*) of 3 and 10 µm, and pore to pore distances (*p* − *p*) of 5 and 10 µm. The other numerical values used for simulations are all listed in Table 2.

**Table 2.** The values of the parameter used for the numerical simulations.

Before presenting the main results, a brief discussion is presented here to investigate the dynamics of NPs injected with an aqueous solution, which substitutes aerosol injection in many relevant experimental studies [8,11,12]. Accordingly, the pore diameter of 10 µm and pore-to-pore distance of 10 µm are taken into consideration, and the deposition rate and transfer rate are compared in Figure 7 for NPs injected via air and media flow into the upper channel of the device. As previously mentioned, the deposition rate is referred to as the proportion of the total released particles that are deposited on the channel substrate. Similarly, the transfer rate defines the ratio between the number of particles transferred to the media channel and the number of released particles. The results obtained demonstrate a noticeable decrease in deposition/transfer rate when media flow is used for the injection of NPs, regardless of their size. To put it clearly, particles follow the fully-developed parabolic profile of media velocity with a lower chance of settling on the substrate or passing through it (data not shown for brevity). Buoyant weight almost fades for submerged particles with a density of the same order as the fluid, which explains the reduction in deposition/transfer rate, especially for large particles. Furthermore, by increasing the fluid viscosity, the drag force dominates the Brownian force, which plays a key role in the deposition of small NPs. Thus, one can draw the conclusion that injecting through an aqueous solution leads to less efficient particle deposition/translocation, and accordingly, air flow injection should be taken into consideration.

#### *4.1. Deposition of Nanoparticles*

This section deals with the sedimentation of airborne NPs in the air channel of the lungon-a-chip device with a focus on the deposition rate and distribution efficiency. Figure 8 depicts the deposition rate against particle diameter for different membrane porosities and air inflow velocities of 0.3 and 1 mm/s. The results show that the deposition rate declines by increasing the particle diameter up to 200 nm regardless of the inlet velocity. The reason behind this decay is related to the dominancy of Brownian diffusion, which is the main particle deposition mechanism for particles of diameter less than 100 nm

according to Stokes–Einstein expression [28,36,37]. On the other hand, as could be deduced from Equations (4)–(6), gravity acceleration becomes more dominant by increasing the particle diameter, while the Brownian effect diminishes. This consequently results in higher deposition rates as the particles enlarge, which can be observed in Figure 8 for diameters larger than 200 nm. A comparison within Figure 8a,b implies an overall decrease in deposition rate by intensifying the inlet flow rate, which makes a larger fraction of NPs leave the device through the air channel outlet. However, both cases follow a similar trend except for particles larger than 700 nm, demonstrating a downtrend of deposition rate for inflow velocity of 0.3 mm/s. By comparing the results of Figure 8a and those of Figure 12a, it is observed that more than 95% of these injected big particles tend to either sediment on the air channel substrate or pass through the perforated membrane with a growing transfer rate against the particle size. Therefore, one can conclude that this unexpected behavior is not linked with blowing NPs off the device but rather with a higher transfer rate into the media channel. Indeed, lower inflow velocity magnitude results in slower axial fluid velocity in the boundary layer above the porous area and prolonged regional residence time, which consequently provide more probability for heavier particles to pass through the membrane. The idea is explicitly illustrated through a series of snapshots in Figure 9, where the trajectory of an initially inert 900 nm particle is compared in the porous region for inflow velocities of 0.3 and 1 mm/s. Such behavior also leads to more dependency of larger NPs' deposition on membrane porosity, which is supported by the divergence of deposition rate curves corresponding with different porosities in Figure 8. , 0.3 mm/s, 1mm/s Media velocity at inlet [23,32] , 0.3 mm/s, 1mm/s Air velocity at inlet [22] Before presenting the main results, a brief discussion is presented here to investigate the dynamics of NPs injected with an aqueous solution, which substitutes aerosol injection in many relevant experimental studies [8,11,12]. Accordingly, the pore diameter of 10 µm and pore-to-pore distance of 10 µm are taken into consideration, and the deposition rate and transfer rate are compared in Figure 7 for NPs injected via air and media flow into the upper channel of the device. As previously mentioned, the deposition rate is referred to as the proportion of the total released particles that are deposited on the channel substrate. Similarly, the transfer rate defines the ratio between the number of particles transferred to the media channel and the number of released particles. The results obtained demonstrate a noticeable decrease in deposition/transfer rate when media flow is used for the injection of NPs, regardless of their size. To put it clearly, particles follow the fullydeveloped parabolic profile of media velocity with a lower chance of settling on the substrate or passing through it (data not shown for brevity). Buoyant weight almost fades for submerged particles with a density of the same order as the fluid, which explains the reduction in deposition/transfer rate, especially for large particles. Furthermore, by increasing the fluid viscosity, the drag force dominates the Brownian force, which plays a key role in the deposition of small NPs. Thus, one can draw the conclusion that injecting through an aqueous solution leads to less efficient particle deposition/translocation, and accordingly, air flow injection should be taken into consideration. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 11 of 18 **Figure 7.** (**a**) Deposition rate and (**b**) transfer rate of the media- and air-induced NPs in the top channel of the lung-on-a-chip device when the fluid velocity is 0.3 mm/s, and both pore diameter and pore-to-pore distance are 10 µm. *4.1. Deposition of Nanoparticles*  This section deals with the sedimentation of airborne NPs in the air channel of the lung-on-a-chip device with a focus on the deposition rate and distribution efficiency. Figure 8 depicts the deposition rate against particle diameter for different membrane porosities and air inflow velocities of 0.3 and 1 mm/s. The results show that the deposition rate declines by increasing the particle diameter up to 200 nm regardless of the inlet velocity. The reason behind this decay is related to the dominancy of Brownian diffusion, which is the main particle deposition mechanism for particles of diameter less than 100 nm according to Stokes–Einstein expression [28,36,37]. On the other hand, as could be deduced from Equations (4)–(6), gravity acceleration becomes more dominant by increasing the particle diameter, while the Brownian effect diminishes. This consequently results in higher deposition rates as the particles enlarge, which can be observed in Figure 8 for diameters

pore to pore distances (−) of 5 and 10 µm. The other numerical values used for simu-

**Parameters Values Descriptions Reference** ℎ 100 µm Height of the channels [32] 2 mm Length of the channels [22] 10 µm, 3 µm Diameter of the membrane pore [32,33] *b* 10 µm Thickness of the membrane [8] ୮ 1180 kg/m3 Density of particles [22] 1000 kg/m3 Density of media [34] 0.718 mPa.s Viscosity of media [34] 1.123 kg/m3 Density of air [35] 0.019 mPa.s Viscosity of air [35]

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 10 of 18

**Table 2.** The values of the parameter used for the numerical simulations.

lations are all listed in Table 2.

**Figure 7.** (**a**) Deposition rate and (**b**) transfer rate of the media- and air-induced NPs in the top channel of the lung-on-a-chip device when the fluid velocity is 0.3 mm/s, and both pore diameter and pore-to-pore distance are 10 µm. for inflow velocities of 0.3 and 1 mm/s. Such behavior also leads to more dependency of larger NPs' deposition on membrane porosity, which is supported by the divergence of deposition rate curves corresponding with different porosities in Figure 8.

**Figure 8.** Deposition rate of airborne NPs of diameter ranging from 10 nm to 900 nm on the substrate of the air channel with different membrane porosities and airflow velocities of (**a**) 0.3 mm/s and (**b**) 1 mm/s. **Figure 8.** Deposition rate of airborne NPs of diameter ranging from 10 nm to 900 nm on the substrate of the air channel with different membrane porosities and airflow velocities of (**a**) 0.3 mm/s and (**b**) 1 mm/s.

**Figure 9.** Trajectory of an initially inert 900 nm particle in the porous region for inflow velocities of (**a**) 0.3 mm/s and (**b**) 1 mm/s. **Figure 9.** Trajectory of an initially inert 900 nm particle in the porous region for inflow velocities of (**a**) 0.3 mm/s and (**b**) 1 mm/s.

In addition to the high deposition rate of the injected solid NPs, their better distribution on the channel substrate is another significant factor in toxicological and pharmaceutical studies to attain an effective biochemical impact on cultured cells. Here, the distribution of NPs is mathematically described in terms of the mean value and standard deviation of their position [38]. In addition to the high deposition rate of the injected solid NPs, their better distribution on the channel substrate is another significant factor in toxicological and pharmaceutical studies to attain an effective biochemical impact on cultured cells. Here, the distribution of NPs is mathematically described in terms of the mean value and standard deviation of their position [38].

Figure 10 depicts dimensionless distribution indexes of deposited NPs on the air channel substrate under two air inflow velocities, which are normalized with respect to channel length here. Since membrane porosity has no significant effect on statistical characteristics of particle distribution due to the regular alignment of pores, only a single case with a pore diameter of 10 µm and pore-to-pore distance of 10 µm is considered. Distribution efficiency can simply be assessed against an ideal uniform deposition with normalized mean value of 0.5 and a standard deviation of 1/√12 (shown with a solid red line in Figure 10) [38]. Figure 10 depicts dimensionless distribution indexes of deposited NPs on the air channel substrate under two air inflow velocities, which are normalized with respect to channel length here. Since membrane porosity has no significant effect on statistical characteristics of particle distribution due to the regular alignment of pores, only a single case with a pore diameter of 10 µm and pore-to-pore distance of 10 µm is considered. Distribution efficiency can simply be assessed against an ideal uniform deposition with normalized mean value of 0.5 and a standard deviation of 1/ √ 12 (shown with a solid red line in Figure 10) [38]. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 13 of 18

**Figure 10.** Normalized (**a**) average position and (**b**) standard deviation of deposited NPs of diameter ranging from 10 nm to 900 on the substrate of the air channel with pore diameter and pore to pore distance of 10 µm under two different fluid velocities of 0.3 mm/s and 1 mm/s. The solid red line corresponds with an ideal deposition. **Figure 10.** Normalized (**a**) average position and (**b**) standard deviation of deposited NPs of diameter ranging from 10 nm to 900 on the substrate of the air channel with pore diameter and pore to pore distance of 10 µm under two different fluid velocities of 0.3 mm/s and 1 mm/s. The solid red line corresponds with an ideal deposition.

Under an airflow velocity of 0.3 mm/s, both small and large particles demonstrate low values of the mean location and standard deviation, indicating a concentrated deposition at the entrance of the channel. Increasing the airflow velocity would not alter the results unless for particles larger than 600 nm or smaller than 50 nm, giving rise to a significantly widened distribution range. Although the uniform distribution achieved by higher inflow velocity

**Figure 11.** Snapshots from transient distribution of NPs with a diameter of 900 nm in the lung-ona-chip device with pore diameter and pore-to-pore distance of 10 µm under the fluid velocity of (**a**)

This section deals with the translocation of airborne NPs into the media channel of the lung-on-a-chip device, mainly analyzing transfer rate and distribution efficiency. Figure 12 compares the transfer rate of NPs for different membrane porosities and inflow velocities. A higher translocation efficiency seems achievable by increasing the former and decreasing the latter. The obtained results show that for constant membrane porosity and inlet velocity, the transfer rate declines by increasing the particle diameter up to 200 nm, whereafter it starts to go up. Higher transfer rates arise from Brownian diffusivity for small particles and dominancy of gravity force, as well as longer residence time in the

0.3 mm/s, and (**b**) 1 mm/s.

*4.2. Translocation of Nanoparticles*

porous regions for large particles.

seems to be desirable, it should be noted that the deposition rate would be unfavorably affected in this case. Therefore, depending on the application, a tradeoff between deposition rate and distribution uniformity should be considered to obtain proper particle size and inflow velocity values. To provide a better insight into the problem, snapshots from the distribution of 900 nm particles are provided in Figure 11. **Figure 10.** Normalized (**a**) average position and (**b**) standard deviation of deposited NPs of diameter ranging from 10 nm to 900 on the substrate of the air channel with pore diameter and pore to pore distance of 10 µm under two different fluid velocities of 0.3 mm/s and 1 mm/s. The solid red line corresponds with an ideal deposition.

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 13 of 18

**Figure 11.** Snapshots from transient distribution of NPs with a diameter of 900 nm in the lung-ona-chip device with pore diameter and pore-to-pore distance of 10 µm under the fluid velocity of (**a**) 0.3 mm/s, and (**b**) 1 mm/s. **Figure 11.** Snapshots from transient distribution of NPs with a diameter of 900 nm in the lung-ona-chip device with pore diameter and pore-to-pore distance of 10 µm under the fluid velocity of (**a**) 0.3 mm/s, and (**b**) 1 mm/s.

#### *4.2. Translocation of Nanoparticles 4.2. Translocation of Nanoparticles*

This section deals with the translocation of airborne NPs into the media channel of the lung-on-a-chip device, mainly analyzing transfer rate and distribution efficiency. Figure 12 compares the transfer rate of NPs for different membrane porosities and inflow velocities. A higher translocation efficiency seems achievable by increasing the former and decreasing the latter. The obtained results show that for constant membrane porosity and inlet velocity, the transfer rate declines by increasing the particle diameter up to 200 nm, whereafter it starts to go up. Higher transfer rates arise from Brownian diffusivity for small particles and dominancy of gravity force, as well as longer residence time in the porous regions for large particles. This section deals with the translocation of airborne NPs into the media channel of the lung-on-a-chip device, mainly analyzing transfer rate and distribution efficiency. Figure 12 compares the transfer rate of NPs for different membrane porosities and inflow velocities. A higher translocation efficiency seems achievable by increasing the former and decreasing the latter. The obtained results show that for constant membrane porosity and inlet velocity, the transfer rate declines by increasing the particle diameter up to 200 nm, whereafter it starts to go up. Higher transfer rates arise from Brownian diffusivity for small particles and dominancy of gravity force, as well as longer residence time in the porous regions for large particles. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 14 of 18

**Figure 12.** Transfer rate of NPs of diameter ranging from 10 nm to 900 to the media channel of the lung-on-a-chip device with different membrane porosities considering equal inflow velocities of (**a**) 0.3 mm/s and (**b**) 1 mm/s for air and media channels. **Figure 12.** Transfer rate of NPs of diameter ranging from 10 nm to 900 to the media channel of the lung-on-a-chip device with different membrane porosities considering equal inflow velocities of (**a**) 0.3 mm/s and (**b**) 1 mm/s for air and media channels.

bution of transferred NPs are assessed at the outlet of media channel utilizing Equations (4) and (5) with representing the number of transferred particles and showing the vertical position of each particle. Here, the average position and standard deviation are normalized with respect to media channel height, and, similar to the previous section, are then assessed against an ideal uniform distribution (shown with the solid red line in Figure 13). Higher average position and smaller standard deviation for large particles in Figure 13 indicate their tendency to accumulate close to the media channel's top surface, where the endothelial cells are normally cultured. On the other hand, smaller particles demonstrate a higher penetration power, which is due to their higher vertical velocity magnitude while crossing the membrane as well as their krelatively strong Brownian force in the medium. Moreover, a more uniform distribution is observed for medium-sized particles with a diameter range of 100–300 nm. For a better vision of the impact of size, one can compare the snapshots from the distribution of 10 nm particles in Figure 14 with those

Distribution of NPs inside the media channel is also of great importance, especially for targeted particle delivery applications. Since submerged particles with a density close

obtained from 900 nm particles in Figure 11a.

Distribution of NPs inside the media channel is also of great importance, especially for targeted particle delivery applications. Since submerged particles with a density close to that of fluid would follow the media velocity profile regardless of inflow velocity, the results are presented only for the inlet velocity of 0.3 mm/s. Accordingly, vertical distribution of transferred NPs are assessed at the outlet of media channel utilizing Equations (4) and (5) with *N* representing the number of transferred particles and *xn* showing the vertical position of each particle. Here, the average position and standard deviation are normalized with respect to media channel height, and, similar to the previous section, are then assessed against an ideal uniform distribution (shown with the solid red line in Figure 13). Higher average position and smaller standard deviation for large particles in Figure 13 indicate their tendency to accumulate close to the media channel's top surface, where the endothelial cells are normally cultured. On the other hand, smaller particles demonstrate a higher penetration power, which is due to their higher vertical velocity magnitude while crossing the membrane as well as their krelatively strong Brownian force in the medium. Moreover, a more uniform distribution is observed for medium-sized particles with a diameter range of 100–300 nm. For a better vision of the impact of size, one can compare the snapshots from the distribution of 10 nm particles in Figure 14 with those obtained from 900 nm particles in Figure 11a. Distribution of NPs inside the media channel is also of great importance, especially for targeted particle delivery applications. Since submerged particles with a density close to that of fluid would follow the media velocity profile regardless of inflow velocity, the results are presented only for the inlet velocity of 0.3 mm/s. Accordingly, vertical distribution of transferred NPs are assessed at the outlet of media channel utilizing Equations (4) and (5) with representing the number of transferred particles and showing the vertical position of each particle. Here, the average position and standard deviation are normalized with respect to media channel height, and, similar to the previous section, are then assessed against an ideal uniform distribution (shown with the solid red line in Figure 13). Higher average position and smaller standard deviation for large particles in Figure 13 indicate their tendency to accumulate close to the media channel's top surface, where the endothelial cells are normally cultured. On the other hand, smaller particles demonstrate a higher penetration power, which is due to their higher vertical velocity magnitude while crossing the membrane as well as their krelatively strong Brownian force in the medium. Moreover, a more uniform distribution is observed for medium-sized particles with a diameter range of 100–300 nm. For a better vision of the impact of size, one can compare the snapshots from the distribution of 10 nm particles in Figure 14 with those obtained from 900 nm particles in Figure 11a.

**Figure 12.** Transfer rate of NPs of diameter ranging from 10 nm to 900 to the media channel of the lung-on-a-chip device with different membrane porosities considering equal inflow velocities of (**a**)

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 14 of 18

0.3 mm/s and (**b**) 1 mm/s for air and media channels.

**Figure 13.** Normalized (**a**) average position and (**b**) standard deviation of transferred NPs of diameters ranging from 10 nm to 900 at the outlet of the media channel in the lung-on-a-chip device with different membrane porosities under the fluid velocity of 0.3 mm/s. The solid red line corresponds with an ideal distribution in the media channel.

In addition to distribution efficiency, the current study investigates the transient evolution of transferred NPs as a significant factor in toxico-pharmacokinetic analyses [39]. First, the relative concentration, which is defined by the ratio of total injected particles that flow out through the outlet of the media channel, is recorded, and then, a numerical differentiation is utilized to obtain the time derivative of the data. After a Gaussian curve fitting, final results are depicted in Figure 15 for different membrane porosities and inlet velocities. Figure 15, which is known as the so-called concentration-time plot [38], provides important details describing the behavior of transferred NPs to the media channel, such as area under each curve, which is equal to the transfer rate (see Figure 12), maximum relative concentration rate (*cmax*), and time to reach this maximal value (*tmax*). For example, it is shown that 10 nm particles demonstrate the highest *cmax*, which is significantly augmented by increasing the inflow velocity. Moreover, it is evident that smaller particles display a faster transmission with a lower *tmax*. Although the overall value of *cmax* grows by increasing the membrane porosity and inflow velocity, it should be noted that this may also considerably change *tmax* depending on the particle size.

with an ideal distribution in the media channel.

**Figure 13.** Normalized (**a**) average position and (**b**) standard deviation of transferred NPs of diameters ranging from 10 nm to 900 at the outlet of the media channel in the lung-on-a-chip device with different membrane porosities under the fluid velocity of 0.3 mm/s. The solid red line corresponds

**Figure 14.** Snapshots from transient distribution of particles with a diameter of 10 nm in the lungon-a-chip device with pore diameter and pore-to-pore distance of 10 µm and fluid velocity of 0.3 mm/s. **Figure 14.** Snapshots from transient distribution of particles with a diameter of 10 nm in the lung-on-a-chip device with pore diameter and pore-to-pore distance of 10 µm and fluid velocity of 0.3 mm/s. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 16 of 18

**Figure 15.** Transient relative concentration rate for NPs with different diameters. (**a**) =3 µm, − = 10 µm; = 0.3 mm/s; (**b**) = 10 µm, −=5 µm, = 0.3 mm/s; (**c**) =3 µm, − = 10 µm, = 1 mm/s; and (**d**) = 10 µm, − = 5 µm, = 1 mm/s. **Figure 15.** Transient relative concentration rate for NPs with different diameters. (**a**) *d* = 3 µm, *p* − *p* = 10 µm; *ua* = 0.3 mm/s; (**b**) *d* = 10 µm, *p* − *p* = 5 µm, *ua* = 0.3 mm/s; (**c**) *d* = 3 µm, *p* − *p* = 10 µm, *ua* = 1 mm/s; and (**d**) *d* = 10 µm, *p* − *p* = 5 µm, *ua* = 1 mm/s.

With a promising potential in recapitulating in vivo conditions, lung-on-a-chip tech-

a sandwiched porous membrane enables particle tracing in air and media phases simultaneously, which is an essential prerequisite for further toxicological and pharmaceutical studies. Despite this remarkable capability, a detailed investigation into the dynamics of NPs in this device has been overlooked by the literature. Here, a comprehensive numerical parametric study was conducted to study the role of various hydrodynamical and geometrical properties on the dynamics of airborne NPs. Accordingly, the impacts of fluid flow velocity, membrane porosity, and particle diameter were examined on the dispersion of the particles in both air and media channels, as well as their deposition and translocation. Consistent with the experimental observation reported in the literature, we numerically showed that the aerosol injection of NPs provided far more efficient deposition/translocation than the aqueous solution. Although very small particles (୮ ൏ 50 nm) as well as very large ones (୮ 700 nm) demonstrated a high deposition rate, they tended to sediment close to the entrance of the channel for low inflow velocities, which could be undesirable for many applications demanding a uniform distribution. On the other hand, a higher inflow velocity provided a more uniform distribution; however, it also resulted in lower deposition and translocation rates. Therefore, it can be concluded that mediumsized particles (500 nm ൏ ୮ ൏ 700 nm) are more appropriate for uniform distribution with lower inflow velocities. It was demonstrated that the translocation of NPs into the media channel became more efficient by reducing the particle size as a factor of Brownian

**5. Conclusions**

#### **5. Conclusions**

With a promising potential in recapitulating in vivo conditions, lung-on-a-chip technology has recently attracted the attention of researchers. This dual-channel device with a sandwiched porous membrane enables particle tracing in air and media phases simultaneously, which is an essential prerequisite for further toxicological and pharmaceutical studies. Despite this remarkable capability, a detailed investigation into the dynamics of NPs in this device has been overlooked by the literature. Here, a comprehensive numerical parametric study was conducted to study the role of various hydrodynamical and geometrical properties on the dynamics of airborne NPs. Accordingly, the impacts of fluid flow velocity, membrane porosity, and particle diameter were examined on the dispersion of the particles in both air and media channels, as well as their deposition and translocation. Consistent with the experimental observation reported in the literature, we numerically showed that the aerosol injection of NPs provided far more efficient deposition/translocation than the aqueous solution. Although very small particles (*d*p < 50 nm) as well as very large ones (*d*p > 700 nm) demonstrated a high deposition rate, they tended to sediment close to the entrance of the channel for low inflow velocities, which could be undesirable for many applications demanding a uniform distribution. On the other hand, a higher inflow velocity provided a more uniform distribution; however, it also resulted in lower deposition and translocation rates. Therefore, it can be concluded that medium-sized particles (500 nm < *d*p < 700 nm) are more appropriate for uniform distribution with lower inflow velocities. It was demonstrated that the translocation of NPs into the media channel became more efficient by reducing the particle size as a factor of Brownian diffusivity and also by increasing the particle residence time in the porous region as a consequence of the gravity force's dominancy in large particles. However, smaller particles showed a better distribution in the media channel than did larger ones. Therefore, the results of this paper propose optimizing particle size as well as inflow velocity and membrane porosity, depending on the objectives of biomedical investigations. Even though the current numerical analysis is carried out with some simplifications, such as neglecting the cultured cell layers, it takes the first step towards a perceptive insight into the dynamics of the NPs in the lung-on-a-chip device. As a reliable complementary method seeking to reduce experimental cost, time and efforts, this study can be beneficial for future preclinical studies on pulmonary pathology associated with inhaled hazardous and toxic environmental particles, as well as therapeutic studies for developing inhalation drug delivery.

**Author Contributions:** Conceptualization, Z.S., P.A. and N.K.; methodology, Z.S.; software, Z.S.; validation, Z.S.; investigation and resources, Z.S.; writing—original draft preparation, Z.S.; writing—review and editing, P.A., N.K. and C.G.; visualization, Z.S.; supervision, P.A., N.K. and C.G.; funding acquisition, Z.S. and P.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.

**Acknowledgments:** Z.S. and P.A. would like to acknowledge the Shahrood University of Technology, which supported this study, and École Polytechnique Fédérale de Lausanne for supporting the experimental part of this study and its specific infrastructure for computational resources. Z.S. and P.A. also want to thank Swiss Government Excellence Scholarship. N.K. acknowledges the fellowship from the Australian Research Council (ARC) Discovery Early Career Research Award (DECRA) DE220100205.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

36. Friedlander, S.K. *Smoke, Dust, Haze: Fundamentals of Aerosol Dynamics*; Oxford University Press: Oxford, UK, 2000.


*3.2. Optimal Infinitesimal Perturbations of Pulsatile Pipe Flow*

We performed a linear non-modal stability analysis of Sexl–Womersley flow at (*Re* = 2400, *Wo* = 8) and amplitudes up to *A* = 1.6, as described in Section 2.3. For *A* ≤ 0.4, the optimal perturbation is the same as for statistically steady pipe flow: an axial two-roll configuration (not shown here). For *A* ≥ 0.5, the optimal perturbation is a streamwise helix and the optimal initial time of perturbation is *t*0 *T* ∈ [0.2, 0.3]. In Figure 3b,c, we show the optimal perturbation for *A* = 1 at the optimal time of perturbation (*t*0) and at the point of maximum energy amplification (*tf*), respectively. Initially, the optimal helical perturbation is localised very close to the pipe wall at the border of the Stokes layer

(δ = 1 √ 2*Wo* ), and it is tilted towards it. Within the rest of the DC phase and the first stages of AC, the perturbation rapidly grows by four orders of magnitude in energy within only 40% of the period. By the time of maximum energy amplification (early AC phase at *t T* = 0.6), the optimal helix has separated from the Stokes layer and moved completely to the outer bulk region, where the stabilising effect of acceleration arrives later. This wall-normal phase lag increases with *Wo* [13] and can be nicely seen for the profiles close around the peak flow rate ( *t T* = 0.5 and 1.0) in Figure 3a. At the end of the process, the helix has been tilted opposite to its original configuration in a process reminiscent of the Orr mechanism. See Xu et al. [20] for more details and for a comprehensive parametric exploration.


**5. Discussion**

The investigation above addresses an emerging issue in urban ES development and planning: the integration of stakeholder perspectives with quantitative estimates of urban ES provision to inform decisions regarding future land cover [4,16,35,44]. The new method that results is then applied to an urban neighborhood through the evaluation of existing public green space and current urban ES supply rates, the collection of resident and city-wide stakeholder priorities, and the translation of these lines of evidence into a set of alternative planting regimes.

#### *5.1. Public Green Space Inventory*

The population density and land cover distribution in the study area is broadly representative of urban residential neighborhoods in the USA [92], with lawn as the dominant vegetated land cover type (Tables S3 and S4), consistent with other urban areas studied in the USA [29] and in Europe [93]. Approximately one-eighth of this area is covered by trees, again consistent with tree coverage of European research sites [93]. At the same time, resident access to public green space is substantially above the USA average [51] (Figure S1), an unusual situation in a neighborhood with below-median income [49,50] that may partly explain the high value residents placed on numerous urban ES (Figure 4). In the evaluation of existing vegetated land cover types, the combined LiDAR and NDVI-based analysis identified tree cover with high accuracy but was less accurate in identifying shrubs (Table S2); it also obscured the co-occurrence of tree cover over herbaceous ground cover and shrubs. We, therefore, recommend the incorporation of further image analysis rules and waveform LiDAR processing in future green cover analyses for the accurate distinction of these land cover types (e.g., [94,95]).

Currently, lawn provides the majority of runoff retention and recreation value in the neighborhood, while trees and woodlands provide the majority of air purification, carbon storage, and cooling, despite their much smaller area. As a result, the question faced by those responsible for increasing urban ES in Eugene (Section 5.2.) is whether existing public lawn in the neighborhood should be replaced, and if so, with what.

#### *5.2. Stakeholder Priorities*

Alignment with local values is known to be critical to the successful planning, development, and management of public urban green space [89,96], and because resident urban ES priorities are locally idiosyncratic (e.g., [40–42,97]), local input is necessary. Accordingly, Eugene's Parks and Open Space Division surveyed thousands of residents and hundreds of city government and operations employees from 2015–2018, culminating in a vision and implementation plan for future green space development [48]. Although the plan focused on the recreational importance of public parks, responding to expressed desires for public gathering space and for additional sports fields [48,90], another prominent goal was to "further the parks and recreation system's capacity to serve as critical infrastructure for clean air, clean water, flood control, carbon sequestration, and climate resilience" [48] (p. 42). Strategies to provide these urban ES have not yet been determined, but an economic assessment of current park value has been completed [98], and a \$50 million tax bond was passed in 2018 to support park operations and development in preparation for future green infrastructure development to expand urban ES delivery [99], showing the timeliness and relevance of the decision-making process we contemplate here.

#### 5.2.1. Resident Priorities

FAN residents valued urban ES highly overall, with most respondents rating 16 of the 17 listed urban ES as "moderately" to "very" important (Table S5). These primarily represented supporting or regulating services, showing the importance of ecological resilience as represented by native plant species, bird and pollinator habitat, stormwater purification, and carbon sequestration [3]. Similarly, the importance of human health and well-being was shown by the priority of air purification. These findings contrasted strongly, however, with results from Paris and Angers, France, and from Porto and Lisbon, Portugal, in which residents valued cultural and provisioning urban ES most highly [41]; Table S5. They also di ffered from the results of six cities in the USA, in which residents viewed native species and pollution mitigation as low priorities [100], and from results of Barcelona, Spain, in which pollination and biodiversity were low priorities as well [97]. Air purification, however, seems to be one of the few urban ES that shows consistently high priority among city dwellers globally, and the high priority of air quality found here is consistent with recent results from China, Portugal, Spain, France, and the Netherlands [39–41,101].

Outdoor recreation was the single cultural urban ES given high priority by FAN residents, consistent with the previous local survey [90] as well as with recent results from Finland and China [42,101]. In contrast, noise reduction and community identity were low priorities among FAN residents, consistent with findings in European cities [41] but unlike those of Guangzhou, China, in which residents ranked noise abatement as a high priority [101]. FAN residents also assigned a low priority to vegetable and fruit production, in contrast with results from cities such as Barcelona [97], revealing an unexpected lack of support given Eugene's promotion of urban farming [102,103]. Still, this result may show that the numerous community and private gardens in the neighborhood have already met this need.

#### 5.2.2. Delphi Analysis

The perspectives of stakeholders who would either be involved in making and influencing green infrastructure development decisions, or responsible for implementing and maintaining any changes, di ffered from each other substantially in the first round but partly converged during the second round, reflecting the exchange of ideas allowed by the Delphi procedure (Supplementary Materials Section S2; Table S8) and showing the promise of this method in reaching consensus among other diverse stakeholder groups. Among the consensus priorities, Delphi participants valued air and stormwater purification most highly, reflecting current urgent needs for these services at the city level (Section 2.1.). Although residential neighborhoods typically contribute fewer total suspended solids and other stormwater pollutants than do roads with heavy tra ffic [104], low-density development typically contributes greater runo ff volume and a ffects more of its watershed than high-density development [105]. Delphi respondents also assigned high value to carbon sequestration, natural beauty, and native plant species (Table 1), but the prospect of dense vegetation raised safety concerns (Table S8), illustrating an internal conflict in the responses that would require eventual resolution.

In contemplating land cover changes, Delphi participants described the recreational and ease-of-maintenance benefits of lawn in neighborhood and city parks; one participant noted that city operations management is currently converting long-established landscape beds into lawn to realize these benefits, adding that lawn is viewed as "neat-looking". These views, combined with the acknowledgment of lawn's minimal delivery of regulating and supporting services, are consistent with those found among municipal land managers in three Swedish cities [106]; however, Delphi participants did not express concern about the costs and greenhouse gas emissions of mowing, as Swedish stakeholders did [106]. Instead, most Delphi participants viewed the "cost and complexity" associated with the installation and maintenance of alternative planting regimes as a greater set of barriers, although they did not elaborate (Table S8). Still, most Delphi participants (>85%) supported the conversion of some public lawn to alternative plantings to increase regulating and supporting services, agreeing that the primary limitation would be the preservation of sightlines along right-of-way planting strips.

#### 5.2.3. Resident and Delphi Comparison

The di fferences between resident and Delphi urban ES priorities appear to have reflected their frames of reference: while Delphi participants contemplated the city scale, residents considered only their own neighborhood, as intended. Illustrating distinctly city-level concerns, one Delphi participant explained:

"What drives public agencies such as our Public Works Department is regulatory requirements—meeting the Federal Clean Water Act and our NPDES Permit requirements (e.g., reducing pollution into local waterways, reducing flooding, and improving air quality). Other items are secondary to the basic welfare and safety of the general public." (Supplementary Materials Section S3, p. 17).

Consistent with the diversity of stakeholder opinion found in related studies, the partial divergence among Delphi and resident priorities supports the inclusion of neighborhood-level input even in city-wide urban ES planning.

A notable discrepancy was found in the stakeholders' view of cost barriers: while Delphi participants agreed that the cost of lawn conversion was an important barrier (Table S7), over 85% of residents expressed willingness to pay for green infrastructure development that would increase urban ES delivery through tax measures or private donations (Figure 5; Table S6). While stated willingness is not a guarantee of future payment, these findings are consistent with survey results regarding willingness to pay for urban ES in Palm Beach, FL, USA, involving biodiversity, outdoor recreation, and flood protection [107], as well as in Wuhan, Changsha, and Nanchang, China, involving climate regulation, cultural services, air quality, erosion prevention, and habitat services [108]. In each of these cases, the willingness to pay exceeded the expected cost per capita despite heterogeneity in the responses. In contrast, residents of Veneto, Italy expressed a willingness to pay for recreational services but felt that biodiversity conservation and landscape quality should be provided without taxes [109], illustrating the importance of personal history and the governmental context in such attitudes [110] and suggesting that urban planners survey their target neighborhoods before assuming particular cost barriers.

#### *5.3. Integration of Stakeholder Priorities with Potential for Local Benefit*

Recent studies have emphasized the need to integrate the ES priorities of multiple stakeholder groups with the quantitative and qualitative potential for the desired ES to be realized [2,16,35,37,44], but none have explored this further. Here, we found that such integration could be readily accomplished by evaluating questions typical of a green space planning process (Section 4.1.) in light of the corresponding stakeholder survey and urban ES supply rate evidence, yielding a set of alternative planting regimes representing each line of evidence. The ability of each vegetated land cover type to provide multiple urban ES then facilitated the creation of an integrated scheme with limited compromise to stakeholder priorities (Section 4.2.). Of central importance, the contrasts among the individual lines of evidence (i.e., quantitative considerations, resident priorities, and Delphi priorities) illustrated the value of consulting all three. In particular, safety and sightline concerns of Delphi participants moderated the extent of woodland required to maximize air and water quality services, while cost concerns limited stormwater planter area; likewise, residents' outdoor recreation priorities limited the conversion of lawn-based playing fields and emphasized widespread inclusion of native flowering shrubs, carrying particular importance because of residents' financial responsibility for any changes as well as their willingness to bear additional costs. As a result, the Integration planting regime is one that could not have been found by any one approach alone, adding strong support to previous recommendations that diverse stakeholder views and quantitative evidence should be considered together in developing urban ES delivery plans [37].