Characterizing viscoelastic materials via ensemble-based data assimilation of bubble collapse observations
Jean-Sebastien Spratt, Mauro Rodriguez, Kevin Schmidmayer, Spencer Bryngelson, Jin Yang, Christian Franck, Tim Colonius
aa r X i v : . [ phy s i c s . c o m p - ph ] A ug Characterizing viscoelastic materials via ensemble-based dataassimilation of bubble collapse observations
Jean-Sebastien Spratt a, ∗ , Mauro Rodriguez a , Kevin Schmidmayer a , Spencer H. Bryngelson a ,Jin Yang b , Christian Franck b , Tim Colonius a a Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA 91125, USA b Department of Mechanical Engineering, University of Wisconsin-Madison, Madison, WI 53706, USA
Abstract
Viscoelastic material properties at high strain rates are needed to model many biological and med-ical systems. Bubble cavitation can induce such strain rates, and the resulting bubble dynamicsare sensitive to the material properties. Thus, in principle, these properties can be inferred viameasurements of the bubble dynamics. Estrada et al. (2018) demonstrated such bubble-dynamichigh-strain-rate rheometry by using least-squares shooting to minimize the difference between sim-ulated and experimental bubble radius histories. We generalize their technique to account foradditional uncertainties in the model, initial conditions, and material properties needed to uniquelysimulate the bubble dynamics. Ensemble-based data assimilation minimizes the computational ex-pense associated with the bubble cavitation model. We test an ensemble Kalman filter (EnKF),an iterative ensemble Kalman smoother (IEnKS), and a hybrid ensemble-based 4D–Var method(En4D–Var) on synthetic data, assessing their estimations of the viscosity and shear modulus of aKelvin–Voigt material. Results show that En4D–Var and IEnKS provide better moduli estimatesthan EnKF. Applying these methods to the experimental data of Estrada et al. (2018) yields simi-lar material property estimates to those they obtained, but provides additional information aboutuncertainties. In particular, the En4D–Var yields lower viscosity estimates for some experiments,and the dynamic estimators reveal a potential mechanism that is unaccounted for in the model,whereby the viscosity is reduced in some cases due to material damage occurring at bubble collapse.
Keywords:
A. dynamics; B. constitutive behaviour; B. viscoelastic material; C. numericalalgorithms; data assimilation
1. Introduction
Measuring the mechanical properties of soft viscoelastic materials at high strain rates (exceed-ing 10 s − ) is a challenging goal of rheometry. These measurements are of particular interest inbiological and medical engineering, where high strain rates occur during impact and blast expo-sure (Bar-Kochba et al., 2016; Sarntinoranont et al., 2012; Meaney and Smith, 2011; Nyein et al.,2010) and during therapeutic ultrasound (Maxwell et al., 2009; Xu et al., 2007; Mancia et al., 2017;Vlaisavljevich et al., 2015; Bailey et al., 2003). Cavitation, which can take place on exposure to ∗ Corresponding author; [email protected]
Preprint submitted to Journal of the Mechanics and Physics of Solids August 12, 2020 ensile waves, induces high strain rates in surrounding materials, and the resulting bubble dynam-ics are sensitive to the adjacent material properties (Estrada et al., 2018). This observation ledEstrada et al. (2018) to propose a high-strain rate rheometer to estimate the viscoelastic proper-ties of polyacrylamide gels through observation of the bubble radius time history during a laser-generated cavitation event in a sample of the material. To infer the viscosity and shear modulus,they developed a least-square fitting technique which minimizes the difference between the measure-ments and the bubble radius history predicted through a mechanical model of a spherical bubblein an assumed neo-Hookean Kelvin–Voigt viscoelastic material.In this paper, we generalize bubble-dynamic rheometry by considering data assimilation (DA)techniques that can potentially improve predictions in uncertainty-prone high-strain-rate regimes.Provided the material properties are observable and the dynamics are sufficiently sensitive to theirvalues, DA provides a solution to this inverse problem that accounts for uncertainty in both themodel and data (Evensen, 2009a,b; Bocquet and Sakov, 2013b; Schillings and Stuart, 2017). Thus,for bubble-dynamics-based rheometry, DA can address additional uncertainties beyond the unknownviscosity and shear modulus, including those associated with measurement noise, additional materialproperties, modeling assumptions, and initial conditions. DA techniques are characterized as filters if the state (and parameters) are updated at each moment based on the prior trajectory or smoothers if the state history and parameters are estimated over a horizon. By considering both filters andsmoothers, we can gain additional insights into whether constant or time-varying parameters bestfit the observed behavior.Three data assimilation methods are considered: an ensemble Kalman filter (EnKF), an iterativeensemble Kalman smoother (IEnKS), and a hybrid ensemble-based 4D–Var method (En4D–Var).EnKF and IEnKS are variations on the classical Kalman Filter (KF) (Kalman, 1960). Their algo-rithms follow the same structure as the KF, which assimilates data at each step of a discrete timeseries. These are Monte Carlo methods that represent the system via an ensemble. The assimilationuses statistics of the ensemble to calculate a sample covariance. This replaces the covariance matrixof the KF and thus the covariance forecast operator, reducing computational cost. En4D–Var is adifferent ensemble-based approach to the assimilation problem. It is a variant of offline 4D varia-tional data assimilation (4D–Var) (Caya et al., 2005; Tr´emolet, 2007), where a guess for the initialcondition is iterated upon to improve the fit to the data over the entire time domain.In section 2, we describe the specific material–bubble-dynamic model, which matches the oneconsidered by Estrada et al. (2018). In section 3, details of the data assimilation algorithms andtheir implementation are provided. We then examine the relative merits of the estimators insection 4 using synthetic data generated by running the model (with additional simulated noise).This allows us to gauge their relative performance for cases with no modeling uncertainties. Next,in section 5, we apply these estimators to established experimental data for polyacrylamide gels.We compare and contrast results with each method and with the the estimates of Estrada et al.(2018). We summarize the conclusions in section 6.
2. Bubble dynamics model
A physical model for the collapsing bubbles is required to characterize the viscoelastic propertiesof surrounding materials. Many spherical bubble dynamics models exist. Of particular relevancehere are those for cavitation in soft materials (Gaudron et al., 2015; Yang and Church, 2005) andspecific numerical methods for solving them (Warnez and Johnsen, 2015; Barajas and Johnsen,2017). We use the model of Estrada et al. (2018), which adopts approximations validated in previ-ous spherical-bubble models (Prosperetti and Lezzi, 1986; Prosperetti et al., 1988; Akhatov et al.,2001; Epstein and Keller, 1972; Keller and Miksis, 1980; Preston et al., 2007). Key assumptionsof this model are that the motion of the bubble and its contents are spherically symmetric, thebubble pressure is spatially uniform (homobaricity), the temperature of the surrounding materialis constant, and that there is no mass transfer of the non-condensible gas across the bubble wall.The Keller–Miksis equation models the radius evolution (Keller and Miksis, 1980), − ˙ Rc ! R ¨ R + 32 − ˙ R c ! ˙ R = 1 ρ Rc ! (cid:18) p b − γR + S − p ∞ (cid:19) + 1 ρ Rc ˙ (cid:18) p b − γR + S (cid:19) , (1)where R is the bubble radius, c the material speed of sound, ρ the material density, p b the bubbleinternal pressure, γ the bubble-wall surface tension, S the stress integral (see (6)), and p ∞ thefar-field pressure. Under the model assumptions, no mass or energy conservation equations areneeded outside the bubble. Furthermore, the conservation of momentum simplifies to an ordinarydifferential equation for the bubble pressure˙ p b = 3 R (cid:20) − κp b ˙ R + ( κ − K ( T ( R )) ∂T∂r (cid:12)(cid:12)(cid:12)(cid:12) r = R + κp b C p,v C p ( k v ( R )) D − k v ( R ) ∂k v ∂r (cid:12)(cid:12)(cid:12)(cid:12) r = R (cid:21) , (2)where κ is the specific heat ratio, K the thermal conductivity, T the gas temperature, C p thespecific heat, D the binary diffusion coefficient, and k v the vapor mass fraction. Subscripts g , v ,and m refer to gas, vapor, and mixture properties. Conservation of energy in the bubble interioryields an equation for the bubble temperature: ρ m C p (cid:18) ∂T∂t + v m ∂T∂r (cid:19) = ˙ p b + 1 r ∂∂r (cid:18) r K ∂T∂r (cid:19) + ρ m ( C p,v − C p,g ) D ∂k v ∂r ∂T∂r , (3)where v m is the radial mixture velocity and ρ m the mixture density. The boundary condition T ( R ) = T ∞ follows from the model assumptions. The radial mixture velocities are computed as v m ( r, t ) = 1 κp b (cid:20) ( κ − K ∂T∂R − r ˙ p b (cid:21) + C p,v − C p,g C p,m D ∂k∂r . (4)Fick’s law describes the mass diffusion process in the bubble. Casting the conservation of massinside the bubble in terms of the mixture density, the vapor mass fraction inside the bubble is ∂k v ∂t + v m ∂k v ∂r = 1 ρ m r ∂∂r (cid:18) r ρ m D ∂k v ∂r (cid:19) . (5)Under the assumption of equilibrium phase change at the bubble wall, the associated boundarycondition at the wall is p v,sat ( T ( R )) = R v k ( R ) ρ m ( R ) T ( R ), where p v,sat is the saturation pressureof the vapor and R v is the gas constant of the vapor.Equations (1), (2), (3), and (5) form a system of equations. This system is evolved in timewith an implicit Runge–Kutta algorithm that uses the trapezoidal rule and backwards differentia-tion at each step (TR–BDF2) (Hosea and Shampine, 1996). The partial differential equations fortemperature and vapor mass fraction are discretized in space via a uniform grid and computedusing second-order-accurate central finite differences. Estrada et al. (2018) showed that the finite-deformation neo-Hookean Kelvin–Voigt model can represent the material response at high strain3ates. In this framework, the material is modeled with a parallel spring (neo-Hookean elastic re-sponse with shear modulus G ) and dashpot (linear viscous response with viscosity µ ). The stressintegral in (1) is S = − G " − (cid:18) R R (cid:19) − R R − µ ˙ RR , (6)where R is the equilibrium bubble radius Gaudron et al. (2015).
3. Data assimilation methods
Two difficulties that drive the choice of data assimilation method are the nonlinearity of thedynamics and large state vector required to discretize the partial differential equations adequately.The former rules out the standard linearized Kalman filter (EKF) (Kalman, 1960) and the latterrenders its direct nonlinear extensions (e.g. the unscented Kalman filter, UKF) computationallyprohibitive. Instead, ensemble-based methods (Evensen, 1994) are considered. They combine com-putational efficiency with nonlinear dynamics by approximating the state covariance via statistics ofa finite (and typically small) ensemble. We consider three specific ensemble methods: an ensembleKalman filter (EnKF), an iterative ensemble Kalman smoother (IEnKS) and a hybrid ensemble-based 4D–Var method.The discretized equations of section 2 are re-written as a nonlinear operator F , and we definethe linear observation function H that maps the state x to measurement space. This yields thediscrete-time dynamical system x k +1 = F ( x k ) + η k , (7) y k = H ( x k ) + ν k , (8)where x k ∈ R d , y k ∈ R n , η k ∼ N (0 , Σ) , ν k ∼ N (0 , Γ) ,F : R d → R d , H : R d → R n . x k is the d -dimensional state comprised of all the dependent variables plus the unknown parameters x = { R, ˙ R, p b , S, T , C , log(Ca) , log(Re) } , (9)which are the bubble-wall radius, velocity, bubble pressure, stress integral, the discretized tem-perature and vapor concentration fields inside the bubble, and the log-Cauchy and log-Reynoldsnumbers, respectively. The Cauchy and Reynolds numbers are defined asCa ≡ p ∞ G and Re ≡ √ ρp ∞ R max µ . (10)These quantities appear in the nondimensionalized model equations of section 2 and the shearmodulus G and viscosity µ can be computed via (10). The forecast operator F maps log(Ca) andlog(Re) to themselves because they are constant in the physical model. Using the logarithm avoids4egative (and thus non-physical) values during the analysis step of the assimilation algorithms(described in sections 3.1, 3.2, and 3.4).The variable y k is the n -dimensional observation (data) at time k . η k is the unknown processnoise (or model error) added to H ( x k ) to retrieve y k . It is assumed to be Gaussian with zero meanand standard deviation Σ. Similarly, ν k is the assumed Gaussian measurement noise added to F ( x k ) to obtain x k +1 , with zero mean and unknown standard deviation Γ. Throughout this study,the only available measurement is the bubble radius. This means that y k is the radius only, andthe observation operator H is the linear map from the state vector to its first element R . In thefollowing, the linear operator H is sometimes represented as the matrix H for clarity ( Hx = H ( x )).The following methods estimate the full state vector x (including parameters of interest log(Ca)and log(Re)) based on observations of y . The EnKF and IEnKS are online (or quasi-online)methods—they optimize the value of x at each time through the simulation. The IEnKS is deemedquasi-online because it uses data from future times as well. The estimation trails the simulationtime by a fixed number of time steps called the lag. Alternatively, the En4D–Var is an offlinemethod, which only optimizes the initial condition for x , taking into account data from the entiretime-domain. The ensemble Kalman filter (Evensen, 1994) represents the probability density function (PDF)for the state of the dynamics through the statistics of an ensemble of q state vectors. It does notrequire an adjoint, or deriving a tangent linear operator to the physical model (Evensen, 2003,2009a). Starting with suitably randomized initial conditions, each ensemble member is propagatedthrough the physical model, and the predictions are then corrected using the ensemble statistics.The ensemble is initialized with a guess for the initial condition x as the mean, and a givencovariance corresponding to the expected error covariance. In practice, each ensemble memberis independently sampled from a normal distribution with mean x and the assumed covariancematrix. Several initialization strategies exist depending on the system and its dynamics. In thepresent case, the nonlinear dynamics render a systematic approach difficult. Instead, the ensemble isinitialized with a covariance corresponding to our best estimate based on simulations, and adjustedthrough trial and error to optimize results. An ensemble size of q = 48 is used for the tests. Thishas shown to give accurate results while keeping computational costs modest. At any given time,the estimated value for the state vector is then taken to be the ensemble average. x k = 1 q q X j =1 x ( j ) k . (11)The filter is broken down into a forecast and an analysis step. In the forecast step, the physicalmodel is used to step the state forward in time with (7). Each representation of the state vector x ( j ) k in the ensemble at time k is propagated through F with ˆ x ( j ) k +1 = F ( x ( j ) k ). Next, in the analysisstep, if an experimental measurement y k +1 is available at the current time step k + 1, then it is usedto correct the forecast. As described in the dynamical system equations, each ensemble memberis mapped to measurement space H ( x k +1 ). The analysis proceeds by minimizing a cost functioninvolving the difference between H ( x ) and the data point y , while accounting for measurementnoise and model error. This cost function marks the key difference between the EnKF and other5nsemble Kalman methods. The EnKF cost function is given by J ( x ) = 12 k y k − H ( x ) k R + 12 k x − ˆ x k k C k (12)= 12 [ y k − H ( x )] T R − [ y k − H ( x )] + 12 [ x − ˆ x k ] T C − k [ x − ˆ x k ] , (13)where R is the measurement noise covariance matrix, which is an input to the algorithm, and C k is the ensemble covariance at time step k . This covariance is defined as C k = A k ( A k ) T , (14)where A k is the state perturbation matrix A k = 1 √ q − h x (1) k − x k , ... , x ( q ) k − x k i . (15)In fact, the minimization does not make use of the covariance matrix directly, but instead uses thestate perturbation matrix and scaled output perturbation matrix HA k defined as HA k = 1 √ q − h y (1) k − y k , ... , y ( q ) k − y k i . (16)The optimization is carried out by finding the minimizer x k satisfying x k = ˆ x k + A k · w k , (17)with w k a correction coefficient. This restricts the solution to the subspace spanned by the scaledperturbation matrix around the prior estimate ˆ x k . The optimization can be restated as w k = argmin w ∈ R q J ( w ) , (18)where J ( w ) = 12 k w k + 12 k y k − H ( ˆ x k ) − HA k ( w ) k R . (19)The solution is unique, and using the Woodbury matrix identity to write the inversion in measure-ment space, can be written as w k = ( HA k ) T [ R + ( HA k )( HA k ) T ] − ( y k − H ( ˆ x k )) . (20)Performing this inversion in the measurement space is in most cases more computationally efficient.Here this is clear, as the measurement space is comprised of only one variable (bubble radius).Once the minimizer is found and the analysis step complete, covariance inflation is applied to theensemble to correct for the (typical) underestimation of the variance with finite (typically small)ensembles (see section 3.3 for details on covariance inflation). Finally, the forecast step can berepeated. 6 .2. The iterative ensemble Kalman smoother Minimizing deviation from data at future times can help to smooth out estimation and fo-cus on longer-term trends. The IEnKS uses information from one or multiple future time stepsin its assimilation, and can thus be an effective tool. While the ensemble initialization andforecast step are the same as that of the EnKF, the difference in the analysis step is twofold.First, the cost function is modified to minimize difference with data at a single or multiple futuretimes (Evensen and van Leeuwen, 2000). The assimilation thus trails the simulation by a numberof time steps (called the lag). Second, it is no longer minimized analytically but iteratively using aGauss–Newton algorithm.The IEnKS method used here is from Bocquet and Sakov (2013a) and Sakov et al. (2012).Bocquet and Sakov (2013b) have shown it to be effective for state and parameter estimation prob-lems with highly nonlinear dynamics. Its cost function can take two forms referred to as ‘singledata assimilation’ (SDA) or ‘multiple data assimilation’ (MDA) (Bocquet and Sakov, 2013a). TheIEnKS–SDA cost function penalizes difference with measurements at a single time step k + L , where L corresponds to the lag of the smoother. It is given by J ( x ) = 12 k y k + L − H ◦ F k → ( k + L ) ( x ) k R + 12 k x − ˆ x k k C k . (21)On the other hand, the IEnKS–MDA cost function minimizes this difference over a data assim-ilation window (DAW) from k + 1 to k + L , and is expressed as J ( x ) = 12 L X i =1 β i k y k + i − H ◦ F k → ( k + i ) ( x ) k R + 12 k x − ˆ x k k C k , (22)where β i are weights attributed to given time steps with P β i = 1. Again, a solution of the form x = ˆ x k + A k · w is sought, but a Gauss–Newton method is used (Bocquet and Sakov, 2013a). Theminimizer w is found by iterating following w ( i +1) = w ( i ) − H − i ) ∆ J ( i ) ( w ( i ) ) , (23)where i is the iteration number, and H is the approximate Hessian H ( j ) = ( q − I + HA T( j ) R − HA ( j ) , (24)where I is the q × q identity matrix. The gradient is given by∆ J ( j ) = − HA T( j ) R − [ y k + L − H ◦ F k + L ← k ( x k )] + ( q − w ( j ) . (25)For the smoother, HA is more complicated than it is for the filter as it involves differenceswith measurements at future time steps. This quantity is akin to a tangent linear operator fromensemble to measurement space and has to be estimated. Following Bocquet and Sakov (2013a), afinite-difference estimate is used: HA ( j ) ≈ α H ◦ F k + L ← k ( x ( j ) k T + α A k ) (cid:18) I − · T q (cid:19) , (26)with scaling factor α ≪ = (1 · · · T a vector of length q . The iteration is repeated untila threshold w ( i +1) − w ( i ) < ǫ , or a fixed number of iterations is reached. Once the optimal value7 opt is obtained, x opt = ˆ x + A k · w opt is calculated and a new ensemble E k is sampled at timestep k with E k = x opt T + p q − A k H − / I . (27)This completes the analysis step. When using the MDA variant, the Hessian and gradient of J arefound with H ( j ) = ( q − I + L X i =1 HA T i β i R − HA i (28)∆ J ( j ) = − L X i =1 HA T i β i R − [ y k + i − H ◦ F k + i ← k ( x i )] + ( q − w ( j ) . (29) While the EnKF and IEnKS may converge, ensemble methods are subject to intrinsic samplingerror (Bocquet, 2011; Luo and Hoteit, 2011). This sampling error results from the finite ensemblesize q used to represent the statistics of a system of often much higher dimension. As van Leeuwen(1999) explains, the EnKF tends to underestimate error variances, particularly for small ensemblesizes. There exist different ways to address this sampling error, but a simple approach is covarianceinflation (Whitaker and Hamill, 2012), where we correct x ( j ) = x + α ( x ( j ) − x ) + λ ( j ) . (30)Here, x denotes the ensemble average, as defined in (11), after the analysis step. Parameters α and λ correspond to multiplicative and additive inflation parameters, respectively.There exist many schemes for multiplicative inflation, the most simple of which is picking ascalar α (usually 1 . ≤ α ≤ . α is found at each time step using α i = 1 + θ (cid:18) σ bi − σ ai σ ai (cid:19) , (31)where σ ai and σ bi are the prior and posterior ensemble standard deviation for the i th element of thestate vector ( α is a vector here), and θ is a scalar (usually 0 . ≤ θ ≤ . α shows, this scheme inflates the covariance more in regions where the analysis led to a large correction.Whitaker and Hamill (2012) test this method and compare it to other approaches, showing thatit performs well. We similarly find that this performs as well or better than a simple scalar α forour tested cases. This RTPS model was used with θ = 0 .
7. Additive covariance inflation was notfound to significantly affect results and introduced some stability issues with larger magnitudes of λ . Therefore, λ = is used. As with ensemble Kalman methods, ensembles can be used with 4D-Var to estimate covarianceempirically, thus reducing computational cost (Gustafsson and Bojarova, 2014; Liu et al., 2008).8he present method (En4D-Var) is a fully offline extension of the IEnKS–MDA method. Again,the ensemble is initialized in the same way as EnKF, but the cost function is here J ( x ) = 12 X k β k k y k − H ◦ F k ← ( x ) k R + 12 k x − ˆ x k k C . (32)The difference with the IEnKS–MDA cost function is the data assimilation window size. Ratherthan minimizing over a few time steps forward and then stepping through time, the minimizationis done over the entire time domain and only the initial state vector is corrected. Each newiteration is initialized with the corrected initial state (including parameters to estimate). Thesame minimization procedure as described in section 3.2 is used. When the minimization hasconverged, a final simulation is run with the forecast model only. In cases where only a fewiterations are necessary, this method reduces computational cost as compared to the IEnKS-MDA.The time dimension is still full included, but each point in time is only assimilated once per iteration.Furthermore, this retains the advantage of ensemble methods. As opposed to classical 4D–Var, thereis no need to linearize the state function and find the tangent linear adjoint operator. This noveladaptation of the IEnKS method is well suited to the present problems given that our interest isthe estimation of material properties which are, at the outset, assumed to be constant.
4. Testing with synthetic data
Synthetic data where the true shear modulus and viscosity are known is generated from themodel (section 2) and used to test the data assimilation methods in a setting where there is nomodeling error. Bubble radius time-history data from the simulation is sampled at 270,000 framesper second to match available experiments. Random Gaussian noise is added to these samplesto mimic experimental data. The standard deviation of this noise is set at σ = 0 .
02, which isgreater than the estimated noise of the experiments. Two polyacrylamide gels were examined withnominal values of shear modulus and viscosity determined by Estrada et al. (2018). For the stiffgel: G stiff = 7 .
69 kPa, µ stiff = 0 .
101 Pa s, and for the soft gel: G soft = 2 .
12 kPa, µ soft = 0 .
118 Pa s.Since similar estimation accuracy was achieved in both cases, we report results for the stiff gel only.The other material properties used are taken from Estrada et al. (2018) and given in table 1. Nouncertainty is added to these parameters in the present study to match their conditions and focuson estimating G and µ . Parameter Value Parameter Value ρ / m c / s p ∞ . γ . × − N / m D . × − m / s κ . C p,g .
62 kJ / kg K C p,v .
00 kJ / kg K A . × − W / m K B . × − W / m K p ref . × kPa T ref T ∞ .
15 K
Table 1:
Model parameters as they follow from Estrada et al. (2018).
9n example simulated radius curve and sampled surrogate measurements (with noise added)with these parameters is shown in figure 1a, plotted against non-dimensional time t ∗ = tR max r p ∞ ρ . (33)With the simulated data, the evolution over time of all variables in the state vector is known. Forexample, bubble-wall velocity, bubble pressure and stress integral are plotted in figure 1b. (a) (b)Figure 1: Simulated bubble radius and noisy sampled data used to test data assimilation methods (a) ,alongside simulated bubble-wall velocity, normalized bubble pressure and stress integral (b) , plotted overnon-dimensional time t ∗ . A set of initial guesses for the shear modulus and viscosity, ranging from 10% to 100% errorfrom the true values, were used to test each method. Table 2 summarizes results for a subset ofthese cases, representing 10, 50 and 100% initial error in G and µ . In each case, ensembles wereinitialized as Gaussian with these erroneous material properties as the mean, and standard devia-tion increasing with increased error. That is, the spread of the initial ensemble was made wider forcases with more error, to account for the increased uncertainty in the initial guess. To match thetests on experimental data in the next section, simulated data is limited to the first three peaks ofthe bubble collapse. This corresponds to approximately 35 points given the initial conditions andframe rate. Estrada et al. (2018) found that limiting the data to this region led to better parameterestimation. Similarly, we find that the model fails to fit the radius measurements after this time.Reasons for this reduced accuracy at later times are discussed in section 5.3.10ethod G estimate (%error) [kPa] µ estimate (%error) [Pa s] Run time [s]Guess 1 8 .
50 (+10%) 0 .
09 (-10%) –EnKF 7 .
234 (5.93%) 0 .
098 (2.61%) 428IEnKS–SDA (lag 1) 7 .
364 (4.24%) 0 .
110 (8.58%) 852IEnKS–MDA (lag 3) 6 .
682 (13.11%) 0 .
100 (0.92%) 8076En4D–Var 7 .
150 (7.03%) 0 .
099 (1.80%) 679Guess 2 3 .
80 (-50%) 0 .
05 (-50%) –EnKF 3 .
988 (48.1%) 0 .
057 (43.1%) 375IEnKS–SDA (lag 1) 4 .
203 (45.4%) 0 .
080 (20.9%) 904IEnKS–MDA (lag 3) 7 .
390 (3.90%) 0 .
086 (15.2%) 9755En4D–Var 7 .
396 (3.82%) 0 .
100 (0.52%) 690Guess 3 15 . .
20 (+100%) –EnKF 13 .
649 (77.5%) 0 .
175 (73.1%) 495IEnKS–SDA (lag 1) 10 .
272 (33.6%) 0 .
142 (40.7%) 800IEnKS–MDA (lag 3) 10 .
078 (31.1%) 0 .
121 (19.9%) 9802En4D–Var 10 .
210 (32.7%) 0 .
118 (16.6%) 611
Table 2:
Comparing accuracy of estimation with 3 different initial guesses for the parameters. Runs wereperformed on a machine with dual 12-core 2.3Ghz processors
Table 2 shows that with a relatively good initial guess with 10% error, the assimilation methodsperform adequately. For example, the EnKF tracks the correct values for shear modulus andviscosity within 6% and 3% respectively. With a moderate initial error of 50%, however, the EnKFlooses accuracy and barely improves on the initial guess. In some cases, the EnKF was observed to beunstable, and the initial ensemble covariance had to be limited to prevent divergence. This limitedthe ability of the filter to estimate the parameters of interest, and thus despite its computationalefficiency, the EnKF is eliminated from further consideration.The IEnKS and En4D–Var performed better than the EnKF for the 50% error case. Theestimation was stable while varying initial conditions and covariance. Still, the lag 1 IEnKS–SDAonly resulted in marginal improvements in the parameter values. The lag 3 IEnKS–MDA, onthe other hand, resulted in further improvement, but at a high computational cost. This cost isassociated with the calculation of the Hessian (see equation (28)) and gradient of the cost function(see equation (29)), which now involves three future time steps. The En4D–Var performs best inthis test case, achieving good estimation with a comparably fast computational time. We note thatwhile the En4D–Var was run for fifteen iterations in each case, the material property estimationconverged by the fifth iterations. Thus, results and run time after five iterations are reported.Estimation results in the case with 50% error are presented in figure 2. Figure 2d shows thesuitability of the En4D–Var: both parameters converge to accurate estimates within a few iterations.Overall, figure 2 also highlights the value of looking over a time horizon. While the EnKF and lag1 IEnKS appear to disbelieve the data too much throughout the run, taking into account multipletimes enables the lag 3 IEnKS–MDA to adjust to new information well, notably around collapse.Indeed, the IEnKS–MDA significantly corrects the viscosity estimate around each collapse, and theshear modulus estimate during the second collapse. Assimilating data from single time-steps appearsto be insufficient given the short time scales of bubble cavitation and limited data. Smoothing overmultiple times far improves performance around collapse points, which, given the IEnKS–MDA11esults, appear to hold the most pertinent data to make the necessary corrections. (a)
EnKF (b)
IEnKS–SDA (lag 1) (c)
IEnKS–MDA (lag 3) (d)
En4D–Var
Figure 2:
Estimation of shear modulus and viscosity with initial guesses of G = 3 . µ = 0 .
05 Pa s(both at 50% error). The estimation is plotted over non-dimensional time t for the EnKF and IEnKSmethods, and over iteration number for the En4D–Var. In the case with 100% error in the initial guess, the relative performances of each method aresimilar to the 50% error case, but the three smoothers stagnate at 20 to 40% errors for µ and G .Weighted by computational expense, the En4D–Var performs best, but the IEnKS–MDA should notbe discarded. Indeed, the time-varying estimation provides additional information about potentiallytime-dependent modeling uncertainties. While the physical model used assumes a constant shearmodulus and viscosity, the quasi-online IEnKS–MDA can uncover potential limitations of thisassumption. This issue is further examined in section 5.3. Ensemble methods carry information about error statistics of the estimated parameters in thefinal ensemble. One way to visualize ensembles is through a histogram, an example of which is shown12 igure 3:
Histogram of the final estimate for log(Ca) with the lag 1 IEnKS and fitted normal curve, wheren is the number of ensemble members at each value of log(Ca). in figure 3 for the logarithm of the Cauchy number with the lag 1 IEnKS–SDA estimator. Despitethe nonlinearity of the model, the tested methods track only the first two statistical moments of anassumed Gaussian filtering or smoothing PDF. Previous works (e.g., (Evensen and van Leeuwen,2000; Yang et al., 2012; Katzfuss et al., 2016)) have discussed that adequate results can still beachieved with a nonlinear model where this assumption must break down to some degree. Ourresults for the IEnKS and En4D–Var results above confirm that this is the case in this example.Figure 4 shows a comparison of the fitted histograms for the methods for the case with 50% initialerror in both parameters. Despite imperfect estimation, the En4D–Var converges significantly morethan other methods given the limited data. The IEnKS–MDA curve displays the least variance ofthe Kalman methods, as expected.Anticipating experimental results, the En4D–Var was run for 10 simulated data sets with thesame ground truth but different (random) noise. Results are shown in figure 5 for shear modulusand viscosity estimates over the data sets. The dashed black lines correspond to the truth, andthe blue line to the mean estimate over the 10 runs. Results across these 10 data sets are fairlyuniform (standard deviation of 0 .
98 kPa for G , 0 .
009 Pa s for µ ), confirming that reliable estimatesare obtained despite noisy measurements across data sets. Figure 6 shows a histogram combiningfinal ensembles for shear modulus to visualize overall results. As each of the 10 ensembles should beapproximately normal, a Gaussian curve is expected when combining them. Figure 6 indeed showsan approximately normal distribution, as does the equivalent histogram for viscosity (as shown insection 5.3 in figure 9a). 13 EnKS-sdaIEnKS-mda EnKFEn4D-VarTruth Initial Guess (a)
EnKFIEnKS-sdaIEnKS-mdaTruth Initial GuessEn4D-Var (b)Figure 4:
Comparing final ensembles for log(Ca) (a) and log(Re) (b) in the case with 50% initial error inboth parameters. (a) (b)Figure 5:
En4D–Var results for G (a) and µ (b) , for ten simulated data sets. igure 6: Histogram for G combining 10 final ensembles for simulated data runs with En4D-Var.
Based on these results with simulated data, given reasonable initial guesses as to the shearmodulus and viscosity, we can confidently expect to estimate both parameters to within 5% usingthe 10 available data sets. Multiple initial guesses can be tested and their fits with experimentalradius histories compared. Therefore, it is straightforward to formulate an initial guess with lessthan 50% error, and thus obtain results comparable to the second test case with simulated data.The En4D–Var is the baseline given its performance. IEnKS–MDA is also tested for its ability toestimate parameters quasi-online.
5. Applied to experimental data
For a more detailed description of the experimental setup for data collection, see Estrada et al.(2018). After the polyacrylamide gel is prepared, each cavitation event is induced with a 6ns pulse ofa “user-adjustable 1–50 mJ, frequency-doubled Q-switched 532 nm Nd:YAG laser”. These cavitationevents are triggered at different locations in the same large batch of polyacrylamide to maximizeuniformity of material properties across experiments. Bubble radius is captured approximatelyevery 3 . µ s, processing 270 000 fps high-speed camera output by subtracting a reference imagefrom each frame and fitting a circle. A few sources of error may be present. Nonuniformity ofthe polyacrylamide gel or discrepancies across data sets could cause the bubble to lose sphericalsymmetry. Laser pulses may also vary slightly across runs, affecting the energy deposited in thesystem and thus initial growth conditions. In practice, a difference in maximum bubble radii wasobserved, with R max = 388 ± µ m across experiments with the stiff gel, and R max = 430 ± µ mwith the soft gel. Ten experimental data sets in the stiff gel were used for the following results, inpart to address this potential lack of uniformity in experimental conditions. G and µ with En4D–Var The noise magnitude in the experimental data is smaller than what was used in the simulateddata with the same data rate. Therefore, if the model is adequate and the noise accurately repre-sented as Gaussian, the IEnKS and En4D–Var should yield comparable or better estimation results15ith the experimental data. As with the simulated data, the assimilation window is limited tothe initial collapse and two subsequent rebounds to match the setup used by Estrada et al. (2018),who estimated in the stiff-gel: G stiff = 7 . ± .
12 kPa and µ stiff = 0 . ± .
023 Pa s. The resultsare compared to theirs. Our estimation is initialized with three different initial guesses, detailed intable 3. Similarly to the surrogate truth data from the last section, initial guesses with 10%, 50%and 100% difference from the Estrada et al. (2018) estimates are chosen.Method G estimate ± σ [kPa] µ estimate ± σ [Pa s] Run time [s]Guess 1 8 .
50 (+10% diff) 0 .
09 (-10% diff) –IEnKS–SDA (lag 1) 7 . ± .
68 0 . ± .
012 2751IEnKS–MDA (lag 3) 7 . ± .
50 0 . ± .
016 9536En4D–Var 7 . ± .
63 0 . ± .
014 609Guess 2 3 .
80 (-50% diff) 0 .
05 (-50% diff) –IEnKS–SDA (lag 1) 4 . ± .
46 0 . ± .
013 2832IEnKS–MDA (lag 3) 6 . ± .
43 0 . ± .
016 10052En4D–Var 6 . ± .
58 0 . ± .
014 585Guess 3 15 . .
20 (+100% diff) –IEnKS–SDA (lag 1) 9 . ± .
76 0 . ± .
014 2871IEnKS–MDA (lag 3) 8 . ± .
52 0 . ± .
015 9222En4D–Var 8 . ± .
58 0 . ± .
016 535
Table 3:
Comparing results of estimation with three different initial guesses for the parameters. Runs wereagain performed on a machine with dual 12-core 2.3Ghz processors
Table 3 summarizes the results from the three different initial material parameters guesses forthe three methods. These estimates correspond to the mean estimate over all 10 experimental datasets. The standard deviation σ of the results is also reported.The En4D–Var estimates for shear modulus and viscosity all fall within the error bounds pro-vided by Estrada et al. (2018). Results are close to theirs in the test cases considered. The estimatesare uniform, with only 8 .
5% difference between viscosity estimates, and 23% difference in the shearmodulus results across the data sets. This larger difference in the shear modulus estimates and inthe associated standard deviations is expected, as we have found the radius curves to be relativelymore sensitive to µ than G . Finally, the average normalized root mean squared error (NRMSE)for bubble radius is low at 2 . × − for guess 1, indicating that a good fit was achieved withthis method (see equation (34) for NRMSE definition). Given that the guess–1 results lead to thesmallest radius error, our initial shear and viscosity modulus estimates are G = 7 . ± .
63 kPa and µ = 0 . ± .
014 Pa s, respectively. An example bubble radius curve is shown in figure 7, for oneof the experimental data sets (data set 10). 16 igure 7:
Radius curve given by En4D–Var estimates and experimental measurements for data set 10.
The standard deviation across the 10 IEnKS–SDA runs was comparable to the En4D–Var.However, the estimates varied significantly based on the initial guess. These ranged from 4 .
32 kPato 9 .
46 kPa for shear modulus, and from 0 .
085 Pa s to 0 .
114 Pa s for viscosity. Except for estimatesfrom guess 1, the shear modulus estimates are also outside the bounds given by Estrada et al. (2018),and the radius fit is significantly worse than that of the En4D–Var, with an average NRMSE of9 . × − .On the other hand, while still worse than that of the En4D–Var, the IEnKS–MDA estimatesare within the Estrada et al. (2018) margin, and the radius fit is better than that of the IEnKS–SDA (NRMSE = 6 . × − ). The IEnKS–MDA thus represents the best tested quasi-onlinemethod, as expected from the simulated data results of section 4. It is important to note here thatthe bubble radius fits were all obtained by re-running simulations with final shear modulus andviscosity estimates and comparing to experimental measurements. This is a fair way to compare theability of each method to estimate these parameters. However, the radius fit obtained online duringassimilation with the IEnKS methods is better (and comparable with the En4D–Var estimates),given that the radius is also being directly corrected at each time-step as part of the state vector.For parameter estimation, the En4D–Var is the best tested method, but the IEnKS–MDA is a goodquasi-online estimator. This is particularly useful for the discussion in section 5.3, where we makeuse of this time-varying estimation. The estimates obtained for the shear modulus and viscosity from the previous section show thatensemble data assimilation methods can be effectively used for estimation of viscoelastic materialproperties. A further look at the results, though, provides more information than simply thisestimate. Examining estimates for each variable across the 10 tested experimental data sets, asshown in figure 8, there appears to be a discrepancy between data sets 3, 4, 5 and the rest for theviscosity. While the shear modulus estimation shows no discernible trend (despite the previouslymentioned larger spread in results), the viscosity data appears to be split between two estimates.The red line in figure 8b shows the mean estimate of data sets 3 to 5 ( µ = 0 .
074 Pa s) and the greenline that of the rest ( µ = 0 .
102 Pa s). 17 a) G estimates (b) µ estimates Figure 8:
En4D–Var estimates for 10 experimental data sets.
Figure 9 compares the histogram obtained when collating the 10 × q final ensemble membersfor viscosity from 10 runs with different simulated data but the same ground truth (Figure 9a) andthe 10 runs done with experimental data (Figure 9b). As discussed in section 4.2, we expect toapproximately retrieve a Gaussian distribution around the estimate, as is the case for the simulatedrun in figure 9a. However, figure 9b shows an apparent bimodal distribution. The lower viscos-ity peak corresponds to the mean estimate of data sets 3 to 5, and the higher peak to that of the rest. (a) Simulated data (b)
Experimental data
Figure 9:
Comparing final combined ensembles for viscosity estimation in simulated and experimental data.
To understand what may be causing this discrepancy of results, it is useful to consider theIEnKS–MDA and its quasi-online estimation of viscosity. Figure 10 shows a comparison between18iscosity estimation for data set 2 (figure 10a) and data set 3 (figure 10b). These are representativeof data sets with a high and low viscosity estimate respectively. They result in estimates of µ =0 .
098 Pa s for data set 2, and µ = 0 .
077 Pa s for data set 3. Comparing the two data sets, itappears that the assimilation begins similarly, correcting to a higher viscosity estimate during thefirst collapse. However, there is a divergence between the behavior of the smoother after eachcollapse point, particularly the second one (around t = 65 µ s). Figure 10a shows a slow decreaseand convergence towards a higher viscosity value, with negligible change at the second collapsepoint. However, this estimate drops sharply after these collapse points in figure 10b. In fact, thedata around each collapse causes the viscosity estimate to sharply drop in data set 3, which doesnot occur in data set 2. This behavior is representative of what is seen in data sets 3 through 5,but does not occur in the rest of the runs. (a) (b)Figure 10: Comparing online estimation of viscosity in data sets 2 (a) and 3 (b) using the IEnKS–MDA.
Given the physical model used, the viscosity should be constant and such drops in the param-eter are not expected. The model alone thus cannot adequately capture the behavior of the gelseen by the IEnKS in these data sets. We can posit that a violent collapse in these data sets iscausing inelastic behavior in the material, and thereupon this perceived change in material prop-erties (Yang et al., 2020). More work will be needed to determine the exact cause, but this couldperhaps result from fracture, damage to the polymer network in the gel, or combustion in the gasphase (Movahed et al., 2016; Kundu and Crosby, 2009; Raayai-Ardakani et al., 2019). Regardlessof physical cause, this time-dependent behavior is not accounted for in the model, but is capturedby the IEnKS–MDA as a drop in the perceived viscosity.Figure 11 compares the normalized root mean squared error (NRMSE) across all data sets, givenby NRMSE = q ( y sim − y exp ) y exp , (34)where y exp is the experimental bubble radius time history, and y sim the simulated time historygiven the estimated material properties (at the corresponding times).19 igure 11: Bar plot of radius normalized root mean squared errors for each data set. Also plotted are theprevious estimate mean NRMSE (mean of all sets), the mean NRMSE for sets 3 to 5, and the final estimatemean NRMSE (mean of all other sets).
Figure 11 shows a higher error in the estimated bubble radius curves fit for data sets 3 through5, which is expected given the heightened model uncertainty in these data sets. Because of thisuncertainty and higher error, we discard these three sets as outliers, which yields the final IEnKS–MDA-informed En4D–Var estimate reported in table 4. Notable are the drop in standard deviationfor viscosity as compared to the previous En4D–Var estimate and the reduced NRMSE.
Estimate G ± σ [kPa] µ ± σ [Pa s] NRMSE Estrada et al. (2018) 7 . ± .
12 0 . ± . . ± .
63 0 . ± .
014 2 . × − Final 7 . ± .
80 0 . ± . . × − Table 4:
Final En4D–Var estimates (discarding three outlier data sets) and standard deviation, along withthe average radius normalized root mean squared error. The previous best estimate corresponds to themean of all 10 data sets, outliers included.
6. Conclusions
Ensemble-based data assimilation was successfully used to estimate the mechanical propertiesof soft viscoelastic materials at high strain rates via observations of bubble collapse. In particular,the ensemble-based 4D–var method (En4D–Var) provided an accurate estimate efficiently, whilethe iterative ensemble Kalman smoother with multiple data assimilation (IEnKS–MDA) reliablyestimated parameters quasi-online. Added benefits of these algorithms include adaptability todifferent numerical or viscoelastic models, and scalability to further parameter estimation, withnegligible computational cost for additional parameters. These methods account for both modeland experimental error, with noisy measurements or inaccuracies in the model having a limitedimpact on estimation. The ability to adjust the entire state vector rather than just the parametersto estimate limits inaccuracies from a poor initial guess. Overall, this represents a viable frameworkfor estimation of mechanical properties of viscoelastic materials.20sing the En4D–Var and IEnKS–MDA together provided information, in this case, in the formof model error. It is hypothesized that the bubble collapses are damaging the polyacrylamidegel in certain test cases, leading to a reduced estimated viscosity after each subsequent collapse,a physical effect not accounted for in the model. Thus, while the En4D–Var provides the bestestimates (especially given its relative computational efficiency), the IEnKS–MDA can provideadditional information about time-dependent modeling errors.
Acknowledgements
This work was supported by the National Institutes of Health [grant number 2P01-DK043881];and the Office of Naval Research [grant numbers N0014-18-1-2625, N0014-17-1-2676].
References
I. Akhatov, O. Lindau, A. Topolnikov, R. Mettin, N. Vakhitova, and W. Lauterborn. Collapse andrebound of a laser-induced cavitation bubble.
Physics of Fluids , 13(10):2805–2819, 2001. doi:10.1063/1.1401810. URL https://doi.org/10.1063/1.1401810 .M. Bailey, V. Khokhlova, O. Sapozhnikov, S. Kargl, and L. Crum. Physical mechanisms of thetherapeutic effect of ultrasound (a review).
Acoust Phys , 49:369–388, 07 2003. doi: 10.1134/1.1591291. URL http://dx.doi.org/10.1134/1.1591291 .E. Bar-Kochba, M. T. Scimone, J. B. Estrada, and C. Franck. Strain and rate-dependent neuronalinjury in a 3d in vitro compression model of traumatic brain injury.
Scientific Reports , 6(1):30550, 2016. doi: 10.1038/srep30550. URL https://doi.org/10.1038/srep30550 .C. Barajas and E. Johnsen. The effects of heat and mass diffusion on freely oscillat-ing bubbles in a viscoelastic, tissue-like medium.
The Journal of the Acoustical Societyof America , 141(2):908–918, Feb 2017. ISSN 0001-4966. doi: 10.1121/1.4976081. URL http://dx.doi.org/10.1121/1.4976081 .M. Bocquet. Ensemble Kalman filtering without the intrinsic need for inflation.
Nonlin-ear Processes in Geophysics , 18(5):735–750, 2011. doi: 10.5194/npg-18-735-2011. URL .M. Bocquet and P. Sakov. An iterative ensemble Kalman smoother.
Quarterly Journal of theRoyal Meteorological Society , 140(682):1521–1535, Oct 2013a. doi: 10.1002/qj.2236. URL http://dx.doi.org/10.1002/qj.2236 .M. Bocquet and P. Sakov. Joint state and parameter estimation with an iterative ensem-ble Kalman smoother.
Nonlinear Processes in Geophysics , 20(5):803–818, Oct 2013b. doi:10.5194/npg-20-803-2013. URL http://dx.doi.org/10.5194/npg-20-803-2013 .A. Caya, J. Sun, and C. Snyder. A comparison between the 4DVAR and the ensemble Kalmanfilter techniques for radar data assimilation.
Monthly Weather Review , 133(11):3081–3094, 2005.doi: 10.1175/MWR3021.1. URL https://doi.org/10.1175/MWR3021.1 .21. Epstein and J. B. Keller. Expansion and contraction of planar, cylindrical, and spherical under-water gas bubbles.
The Journal of the Acoustical Society of America , 52(3B):975–980, 1972. doi:10.1121/1.1913203. URL https://doi.org/10.1121/1.1913203 .J. B. Estrada, C. Barajas, D. L. Henann, E. Johnsen, and C. Franck. High strain-rate soft material characterization via inertial cavitation.
Journal of the Mechanics andPhysics of Solids , 112:291–317, 2018. doi: https://doi.org/10.1016/j.jmps.2017.12.006. URL .G. Evensen. Sequential data assimilation with a nonlinear quasi-geostrophic modelusing Monte Carlo methods to forecast error statistics.
Journal of Geophysi-cal Research: Oceans , 99(C5):10143–10162, 1994. doi: 10.1029/94JC00572. URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94JC00572 .G. Evensen. The ensemble Kalman filter: Theoretical formulation and practical implementa-tion.
Ocean Dynamics , 53(4):343–367, Nov 2003. doi: 10.1007/s10236-003-0036-9. URL http://dx.doi.org/10.1007/s10236-003-0036-9 .G. Evensen.
Data Assimilation: The Ensemble Kalman Filter . Springer, 2009a.G. Evensen. The ensemble Kalman filter for combined state and parameter estimation.
IEEEControl Systems Magazine , 29(3):83–104, June 2009b. doi: 10.1109/MCS.2009.932223. URL http://dx.doi.org/10.1109/MCS.2009.932223 .G. Evensen and P. J. van Leeuwen. An ensemble Kalman smoother for nonlinear dynamics.
MonthlyWeather Review , 128:1852–1867, Jun 2000.R. Gaudron, M. T. Warnez, and E. Johnsen. Bubble dynamics in a viscoelastic medium withnonlinear elasticity.
Journal of Fluid Mechanics , 766:54–75, Jan 2015. doi: 10.1017/jfm.2015.7.URL http://dx.doi.org/10.1017/jfm.2015.7 .N. Gustafsson and J. Bojarova. Four-dimensional ensemble variational (4d-en-var) data assimilationfor the high resolution limited area model (hirlam).
Nonlinear processes in geophysics , 21(4):745–762, 2014.M. Hosea and L. Shampine. Analysis and implementation of TR–BDF2.
Applied NumericalMathematics , 20(1):21–37, 1996. doi: https://doi.org/10.1016/0168-9274(95)00115-8. URL .R. E. Kalman. A new approach to linear filtering and prediction problems.
Journal of BasicEngineering , 82(1):35–45, 1960.M. Katzfuss, J. R. Stroud, and C. K. Wikle. Understanding the ensemble Kalman filter.
The American Statistician , 70(4):350–357, 2016. doi: 10.1080/00031305.2016.1141709. URL https://doi.org/10.1080/00031305.2016.1141709 .J. B. Keller and M. Miksis. Bubble oscillations of large amplitude.
The Journal of the AcousticalSociety of America , 68(2):628–633, Aug 1980. ISSN 0001-4966. doi: 10.1121/1.384720. URL http://dx.doi.org/10.1121/1.384720 . 22. Kundu and A. J. Crosby. Cavitation and fracture behavior of polyacrylamide hydrogels.
Soft Mat-ter , 5:3963–3968, 2009. doi: 10.1039/B909237D. URL http://dx.doi.org/10.1039/B909237D .C. Liu, Q. Xiao, and B. Wang. An Ensemble-Based Four-Dimensional Variational Data AssimilationScheme. Part I: Technical Formulation and Preliminary Test.
Monthly Weather Review , 136(9):3363–3373, 09 2008.X. Luo and I. Hoteit. Robust ensemble filtering and its relation to covariance inflation in theensemble Kalman filter.
Monthly Weather Review , 139(12):3938–3953, 2011. doi: 10.1175/MWR-D-10-05068.1. URL https://doi.org/10.1175/MWR-D-10-05068.1 .L. Mancia, E. Vlaisavljevich, Z. Xu, and E. Johnsen. Predicting tissue susceptibility to me-chanical cavitation damage in therapeutic ultrasound.
Ultrasound in Medicine & Biology ,43(7):1421 – 1440, 2017. doi: https://doi.org/10.1016/j.ultrasmedbio.2017.02.020. URL .A. D. Maxwell, C. A. Cain, A. P. Duryea, L. Yuan, H. S. Gurm, and Z. Xu. Noninvasive throm-bolysis using pulsed ultrasound cavitation therapy – histotripsy.
Ultrasound in Medicine & Bi-ology , 35(12):1982 – 1994, 2009. doi: https://doi.org/10.1016/j.ultrasmedbio.2009.07.001. URL .D. F. Meaney and D. H. Smith. Biomechanics of concussion.
Clinics in SportsMedicine , 30(1):19 – 31, 2011. doi: https://doi.org/10.1016/j.csm.2010.08.009. URL . Concussion inSports.P. Movahed, W. Kreider, A. D. Maxwell, S. B. Hutchens, and J. B. Freund. Cavitation-induced dam-age of soft materials by focused ultrasound bursts: A fracture-based bubble dynamics model.
TheJournal of the Acoustical Society of America , 140(2):1374–1386, 2016. doi: 10.1121/1.4961364.URL https://doi.org/10.1121/1.4961364 .M. K. Nyein, A. M. Jason, L. Yu, C. M. Pita, J. D. Joannopoulos, D. F. Moore, and R. A.Radovitzky. In silico investigation of intracranial blast mitigation with relevance to militarytraumatic brain injury.
Proceedings of the National Academy of Sciences , 107(48):20703–20708,2010. doi: 10.1073/pnas.1014786107. URL .A. T. Preston, T. Colonius, and C. E. Brennen. A reduced-order model of diffusive effects on thedynamics of bubbles.
Physics of Fluids , 19(12):123302, 2007. doi: 10.1063/1.2825018. URL https://doi.org/10.1063/1.2825018 .A. Prosperetti and A. Lezzi. Bubble dynamics in a compressible liquid. part 1. first-order the-ory.
Journal of Fluid Mechanics , 168:457–478, 1986. doi: 10.1017/S0022112086000460. URL http://dx.doi.org/10.1017/S0022112086000460 .A. Prosperetti, L. A. Crum, and K. W. Commander. Nonlinear bubble dynamics.
The Jour-nal of the Acoustical Society of America , 83(2):502–514, 1988. doi: 10.1121/1.396145. URL https://doi.org/10.1121/1.396145 .S. Raayai-Ardakani, D. R. Earl, and T. Cohen. The intimate relationship between cavita-tion and fracture.
Soft Matter , 15:4999–5005, 2019. doi: 10.1039/C9SM00570F. URL http://dx.doi.org/10.1039/C9SM00570F . 23. Sakov, D. S. Oliver, and L. Bertino. An iterative EnKF for strongly nonlinear systems.
Monthly Weather Review , 140(6):1988–2004, Jun 2012. doi: 10.1175/mwr-d-11-00176.1. URL http://dx.doi.org/10.1175/MWR-D-11-00176.1 .M. Sarntinoranont, S. J. Lee, Y. Hong, M. A. King, G. Subhash, J. Kwon, and D. F.Moore. High-strain-rate brain injury model using submerged acute rat brain tissueslices.
Journal of Neurotrauma , 29(2):418–429, 2012. doi: 10.1089/neu.2011.1772. URL https://doi.org/10.1089/neu.2011.1772 . PMID: 21970544.C. Schillings and A. M. Stuart. Analysis of the ensemble Kalman filter for inverse problems.
SIAMJournal on Numerical Analysis , 55(3):1264–1290, 2017.Y. Tr´emolet. Incremental 4D–Var convergence study.
Tellus A: Dynamic Meteorologyand Oceanography , 59(5):706–718, 2007. doi: 10.1111/j.1600-0870.2007.00271.x. URL https://doi.org/10.1111/j.1600-0870.2007.00271.x .P. J. van Leeuwen. Comment on “Data assimilation using an ensem-ble Kalman filter technique”.
Monthly Weather Review , 127(6):1374–1377,1999. doi: 10.1175/1520-0493(1999)127 h i https://doi.org/10.1175/1520-0493(1999)127<1374:CODAUA>2.0.CO;2 .E. Vlaisavljevich, K.-W. Lin, M. T. Warnez, R. Singh, L. Mancia, A. J. Putnam, E. Johnsen,C. Cain, and Z. Xu. Effects of tissue stiffness, ultrasound frequency, and pressure on histotripsy-induced cavitation bubble behavior. Physics in medicine and biology , 60(6):2271–2292, 2015. doi:10.1088/0031-9155/60/6/2271. URL http://dx.doi.org/10.1088/0031-9155/60/6/2271 .M. T. Warnez and E. Johnsen. Numerical modeling of bubble dynamics in viscoelastic mediawith relaxation.
Physics of Fluids , 27(6):063103, Jun 2015. doi: 10.1063/1.4922598. URL http://dx.doi.org/10.1063/1.4922598 .J. S. Whitaker and T. M. Hamill. Evaluating methods to account for system errors in ensemble dataassimilation.
Monthly Weather Review , 140(9):3078–3089, 2012. doi: 10.1175/MWR-D-11-00276.1. URL https://doi.org/10.1175/MWR-D-11-00276.1 .Z. Xu, M. Raghavan, T. L. Hall, C. Chang, M. Mycek, J. B. Fowlkes, and C. A. Cain. High speedimaging of bubble clouds generated in pulsed ultrasound cavitational therapy - histotripsy.
IEEETransactions on Ultrasonics, Ferroelectrics, and Frequency Control , 54(10):2091–2101, 2007.J. Yang, H. C. Cramer, and C. Franck. Extracting non-linear viscoelastic material propertiesfrom violently-collapsing cavitation bubbles.
Extreme Mechanics Letters , 39:100839, 2020. doi:10.1016/j.eml.2020.100839. URL http://doi.org/10.1016/j.eml.2020.100839 .S.-C. Yang, E. Kalnay, and B. Hunt. Handling nonlinearity in an ensemble Kalman filter: Experi-ments with the tree-variable Lorenz model.
Montly Weather Review , 140(8):2628–2646, 2012.X. Yang and C. C. Church. A model for the dynamics of gas bubbles in soft tissue.
The Journalof the Acoustical Society of America , 118(6):3595–3606, Dec 2005. doi: 10.1121/1.2118307. URL http://dx.doi.org/10.1121/1.2118307http://dx.doi.org/10.1121/1.2118307