Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Continuum modeling techniques and their application to the physics of soil liquefaction and dissipative vibrations
(USC Thesis Other)
Continuum modeling techniques and their application to the physics of soil liquefaction and dissipative vibrations
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Continuum Modeling T echniques and
Their Application t o t he Ph ysics of Soil
Liquef action and Dissipativ e V ibr ations
Daniel Lak eland
U niv ersity of Sout her n Calif or nia
December 2013
D edication
This dissertation is dedicated t o m y dear wif e Dr . F r ancesca Mariani, who
supported me during m y g r aduate s tudies and pro vided t he v as t ma jority of
funding f or t his research. She has a deep and abiding lo v e of scientific dis-
co v er y , a s trong sense of what t he tr ul y im portant ques tions are in her field,
and is extremel y g enerous wit h her expertise. Ma y y ou continue t o ha v e
an endless suppl y of g r aduate, under g r aduate, and v olunteer researchers t o
enrich.
Thank y ou t o m y tw o doodlers Nicolas and Alexander , t he mos t im por -
tant products of m y research.
I w ould lik e t o t hank Dr . R og er Ghanem and Dr . Am y R echenmacher
who advised me during m y g r aduate career , and especiall y go t me s tarted
on t he right f oo t, and helped me t o f ocus on t he essential ques tions t hat led
t o m y bigg es t disco v eries.
T o Dr . P aul N e wt on, who seemed t o alw a ys unders tand where I w as
going and appreciated t he nature of m y w or k. His email sa ying “I lik e t he
w a y y ou are approaching it from t he g round up. k eep w or king t hat w a y
- sometimes it’ s har d t o do, and sometimes under appreciated, but I’m a
supporter!” pro vided me wit h enough encour ag ement f or o v er a y ear of
g rinding out results.
T o Dr . Y ehuda Ben-Zion whose research g roup in t he Eart h Sciences
s tudies eart hquak es bo t h in deep details and across t he whole breadt h of
related phenomena. Thank y ou f or y our encour ag ement and man y ex cel-
lent discussions on wide r anging t opics and f or t he en vironment y ou create
of bo t h nurturing and challenging s tudents.
T o Dr . Andre w Gelman whose blog and book s introduced me t o t he
pr actical application of Ba y esian S tatis tics t o scientific data anal ysis, a tech-
nique especiall y suited t o s tatis tical anal ysis in t he context of com plex ph ys-
ical models and measurement errors. Also, f or t he Zombies.
ii
T o Dr . Edw ar d N elson whose IS T sys tem f or nons tandar d anal ysis is suf-
ficientl y sim ple t hat it can be used f or pr actical pur poses b y mere mortals.
T o Dr . Dennis Hazelett who I’ v e kno wn since 5t h g r ade, who pro vided
me wit h at leas t f our or fiv e dis tr acting projects in Biology and Bioinf or matics
during m y PhD, which go t me m y firs t peer re vie w ed publications and g a v e
me f ocused projects where I could refine techniques and sof tw are kno w l-
edg e t hat ultimatel y contributed t o m y ability t o finish m y o wn research.
Hopefull y none of m y g r aphs need t o be tur ned upside do wn.
T o m y mo t her , who alw a ys supported m y ner dy tendencies, and some-
ho w manag ed t o affor d com puters f or me while I w as g ro wing up, e v en
when t he y cos t six mont hs of discretionar y spending, and m y f at her , who
f or all his troubles in lif e w as alw a ys proud of m y accom plishments, I’m
happ y t hat t his w or k v alidates t he im portance of his contribution t o t he UC
Da vis NEES Geo technical Centrifug e project.
Also, I w ould lik e t o t hank Dr . Guy Miller , a colleague at Barr a Inc. in t he
late 1990’ s who firs t go t me t hinking about t he beauty of ph ysics and t he jo y
of modeling and data anal ysis, and Dr . Michael Bishop whose under g r adu-
ate classes in t he Philosoph y Of Science g a v e me an appreciation f or g etting
right t o t he heart of t he matter and no t f ooling around on t he peripher y .
F inall y , t his research w ould ha v e been com putationall y painful if it w ere
no t f or t he w or k of man y contribut ors t o t he F ree Sof tw are projects R, Max-
ima, and J A GS and t he contribut ors t o t he Maxima mailing lis t who pro vide
insight int o an under -appreciated piece of free sof tw are.
iii
Contents
Dedication ii
Chapter 1 Pur pose and Outline 1
Chapter 2 Backg round 3
Chapter 3 T echniq ues and Philosophical Considerations f or Con-
tinuum Models 8
3.1 Continuum Model as S tatis tical Mechanics . . . . . . . . . . . 9
3.2 N on-dimensionalization, Multiple Scales, IS T and N ons tan-
dar d N umbers . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3 The Model Cons tr uction T em plate: Measurement and Mo lecules 17
3.4 The s tatis tical nature of t he model, f alsifiability , par ameter es-
timation, and goodness of fit . . . . . . . . . . . . . . . . . . . 23
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Chapter 4 The Deg ree of Drainag e Dur ing Liq uef action 26
4.1 Back g round . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2 Model Cons tr uction . . . . . . . . . . . . . . . . . . . . . . . . 31
4.3 Qualitativ e R esults . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.4 The role of g r ain sizes and dr ag f orces f or liquef action conse-
quences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.5 The problems wit h small scale em pirical tes ts . . . . . . . . . 56
4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Chapter 5 Dissipation of W a v es in a Molecular T est Problem 61
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.2 The Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
iv
Contents
5.3 The Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.4 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.5 Equiv alent Continuum Model . . . . . . . . . . . . . . . . . . 72
5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
Chapter 6 Com par ing Dissipating W a v e Model to Det ailed T ime-
ser ies 76
6.1 Issues W it h F itting Ba y esian Models t o Dynamics . . . . . . . 81
Chapter 7 Future Directions f or Continuum Modeling of Sand Grains
Dur ing Liq uef action 86
7.1 The Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.2 A sliding t hought experiment . . . . . . . . . . . . . . . . . . 90
7.3 Basic Mechanical Models Com pared . . . . . . . . . . . . . . . 93
7.4 The R ole of Entrop y -lik e Quantities in Modeling Sand . . . . . 96
Chapter 8 Conclusion 99
Bibliog raph y 101
v
Chapter 1
Purpose and Outline
I go among t he fields and catc h a g lim pse of a s t oat or a fieldmouse peeping
out of t he wit her ed gr ass — The cr eatur e hat h a pur pose and its e y es ar e
bright wit h it — I go amongs t t he buildings of a city and I see a man
hurr ying along—t o what?
— John K eats
The research in t his dissertation encom passes f our main t hr us ts: mat hemat-
ical modeling techniques f or building continuum models of sys tems t hat are
explicitl y ackno w ledg ed t o be non-continuous at t he fines t scale o f obser v a-
tion, application of a continuum model of w ater flo w in porous media t o t he
problem of “dr ained v s undr ained” processes in soil liquef action, an exam-
ple model f or t he dissipation of w a v es in a sys tem of molecules decoupled
from an y r adiativ e process, where t he t o tal ener gy s ta ys cons tant, and t he ap-
plication of Ba y esian S tatis tics t o identifying t he par ameters of t hese ph ysical
models.
All of t his w or k arises in t he context of an original project t o anal yze ho w
soil liquef action occurs, and t o incor por ate ph ysical principles such as en-
er gy con v ersion and entrop y g ener ation int o t he ph ysical descrip tion.
The w or k on soil liquef action has been submitted t o Proceedings of t he
R o y al Society A (Lak eland, R echenmacher , and Ghanem n.d. ). A paper is
being prepared on unifying a spectr um of continuum modeling techniques
under t he umbrella of anal ysis via Inter nal Set Theor y (IS T), and t he model
f or dissipativ e w a v es is also being prepared f or submission.
The dissertation proceeds b y firs t describing t he initial mo tiv ation re-
lated t o soil liquef action, and t he back g round kno w ledg e t hat led t o t he re-
search prog r am. Then, I discuss a g ener alized modeling technique using t he
1
mat hematics of IS T which is applicable t o building and inter preting contin-
uum models of explicitl y non-continuous phenomena. F ollo wing t hat back -
g round, a continuum model is presented f or t he fluid flo w t hrough porous
media arising during eart hquak e liquef action and based on what might be
called classical continuum modeling in v ol ving Darcy’ s la w , wit hout explicit
ref erence t o t he IS T based technique. F ollo wing t hat, a model is de v eloped
using t he IS T technique explicitl y which ultimatel y leads in a s tr aightf or -
w ar d w a y t o a dynamic model of dissipativ e w a v es in molecular scale bars,
and a com parison of t he continuum type model t o t he s tatis tics of a detailed
Molecular Dynamics (MD) simulation. F inall y t he com parison t o MD is per -
f or med b y Ba y esian anal ysis t o deter mine t he im plications of t he data f or t he
unkno wn par ameters.
2
Chapter 2
Background
I r emember when our whole island w as shaken wit h an eart hq uake some
y ear s ago, t her e w as an im pudent mount ebank who sold pills whic h (as
he t old t he countr y people) w er e v er y good ag ains t an eart hq uake.
— Joseph A ddison The T atler
Alt hough t he research encom passed in t his dissertation r ang es o v er a wide
sw at h of scales, and includes mat hematical, s tatis tical, mechanical, and philo-
sophical t opics, all of t he research w as ultimatel y mo tiv ated b y initial in-
quiries int o t he w a y in which soil liquefies during an eart hquak e. W it h
apologies f or t he pun, t his t opic has been considered as “settled” in t he
Geo technical field t o t he extent t hat textbook definitions and descrip tions
are consis tent and ha v e been f or man y y ears. The onl y problem is t hat t hese
textbook descrip tions do no t offer much explanat or y po w er in ter ms of w ell
kno wn ph ysical principles, especiall y t he t her modynamics of fluids t hat de-
scribes ho w fluid pressure is g ener ated. The back g round giv en here mo ti-
v ates and inf or ms t he in-dep t h ph ysical anal ysis giv en in future chap ters.
Liquef action of sandy soils is a serious risk in eart hquak e prone regions.
Soil liquef action is a r ang e of phenomena all related t o t he inter action be-
tw een fluid in t he pores of a g r anular soil and t he mo tion of t he g r ains.
The w ors t consequences occur when g r ains cease t o maintain tight contact
and t he mixture is able t o flo w lik e a slurr y , resulting in lar g e def or mations
of t he soil. During m y research t he f oo tag e tak en b y Brent K ooi of lique-
f action in a par k during t he T ohoku-Oki eart hquak e in Japan pro vided a
s triking real-time exam ple of ho w quickl y eart hquak e liquef action can occur
( http://www.youtube.com/watch?v=rn 3 oAvmZY 8 k ). This video pro vides a
clear timescale f or t he liquef action process as t he g round floats on a la y er of
3
w ater and t he w ater rises all t he w a y t o t he surf ace wit hin t he firs t 40 sec-
onds of an eart hquak e whose shaking las ted se v er al minutes in man y places
(Nyquis t 2012 ).
So f ar models of liquef action in engineering pr actice ha v e in v ol v ed pri-
maril y s tatis tical prediction r ules in v ol ving reg ression or cur v e fitting t o case
s tudies and labor at or y tes ts (cf. (Idriss and Boulang er 2006 )). Dynamic mod-
els of liquef action ha v e been built based on continuum models wit h s tress-
s tr ain relationships defined based on phenomenological h ys teresis models
in v ol ving “backbone cur v es” (cf. (Liang 1995 )). Se v er al in v es tig at ors ha v e
sho wn t hat w ater pressure build up is highl y correlated wit h ener gy dis-
sipation, as originall y proposed b y N emat-N asser and Shok ooh ( 1979 ). In
particular , labor at or y and field obser v ations inf or m t he w or k b y F igueroa
et al. ( 1994 ) and Da vis and Berrill ( 2001 ).
Alt hough t he com plexity of t he p henomenon is widel y ackno w ledg ed,
so f ar , t he ph ysical unders tanding of soil liquef action has relied on a core set
of assum p tions and unders tanding t hat is consis tent across t he field. A crit-
ical concep t is t he e ffectiv e s tress, t he t o tal s tress across a section minus t he
portion pro vided b y w ater pressure. An effectiv e s tress of zero means t hat
g r ain-g r ain contacts ha v e disappeared and t his is more or less t he definition
of a full y liquefied soil.
T extbook Definitions The process of liquef action has a consis tent textbook
descrip tion in v ol ving an assum p tion of “undr ained” pore pressure increase.
F or exam ple from Holtz and K o v acs “ The loose sand trieds [sic] t o densify
during shear and t his tends t o squeeze t he w ater out of t he pores. N or mall y ,
under s tatic loading, t he sand has sufficient per meability so t he w ater can
escape and an y induced pore w ater pressures can dissipate. But in t his situ-
ation because t he loading occurs in such a short time, t he w ater doesn ’ t ha v e
time t o escape and t he pore w ater pressure increases. ” (Holtz and K o v acs
1981 , p. 245)
Alter nativ el y from Kr amer : “ The tendency f or dr y cohesionless soils t o
densify under bo t h s tatic and cy clic loading is w ell kno wn. When cohension-
less [sic] soils are satur ated, ho w e v er , r apid loading occurs under undr ained
conditions, so t he tendency f or densification causes ex cess pore pressures t o
increase and effectiv e s tresses t o decrease. ” (Kr amer 1996 , p. 349).
Bo t h descrip tions in v ol v e an ana logy betw een t he dr y sand, which den-
sifies during shearing, and t he satur ated sand which is assumed t o be iso-
4
v olumetric and ins tead increases in w ater pressure. Holtz and K o v acs offer
a picture of a tablet op tank where a quick blo w t o t he side of t he tank in-
s tantl y liquefies t he sand, and piezometers sho w significant w ater pressure
increase o v er s tatic conditions. (Holtz and K o v acs 1981 , p. 243 fig 7.12).
U nf ortunatel y t hese s tandar d descrip tions do no t offer a clear ph ysical
picture of what occurs t o cause t he pressure increase, but t he y bo t h ag ree
on t he undr ained nature of t he e v ent. The tablet op demons tr ation of Holtz
and K o v acs ho w e v er occurs wit h im pact loading las ting per haps 10 or 100
milliseconds, whereas eart hquak e loading occurs o v er a timescale of 10 t o
100 seconds or more. The maximum w ater flux which can occur on a mil-
lisecond time scale is significantl y different from t he mass flux possible o v er
a 10 second time scale.
Problems Wit h The T extbook Descr iption Alt hough t he s tandar d descr ip-
tion is quite sugg es tiv e, t here are man y open ques tions. Ho w are labor at or y
conditions and field conditions similar or different? Does pore w ater pres-
sure build-up cause loss of effectiv e s tress, or does loss of effectiv e s tress
cause pore w ater pressure build up? What is t he role of v olume chang es
and mass flo w? K okusho and K ojima offer an exam ple labor at or y experi-
ment where w ater mig r ation is clear l y im portant and w ater film f or mation is
t he primar y explanation f or soil s trengt h f ailure (K okusho and K ojima 2002 ;
K okusho 1999 ). These experiments in v ol v e a column of soil in t he vicinity of
1 t o 2 meters high carefull y la y ered, and ex cited b y a blo w t o t he bo tt om of
t he co lumn. Af ter a brief tr ansient period, under certain circums tances t he
bo tt om la y ers of soil settle do wn w ar d and a long-las ting w ater film f or ms
betw een t he bo tt om la y er and upper la y ers. Clear l y , w ater has mig r ated up-
w ar ds while soil has mig r ated do wn w ar ds, all wit hin a timescale of a sec-
ond or so. Such beha vior is inconsis tent wit h t he concep t of an “undr ained”
e v ent. Ho w can w e explain t hese phenomena in a ph ysicall y consis tent w a y?
The W a y F or w ard The insight t hat w a v e ener gy dissipation w as associ-
ated wit h soil w ater pressure chang es as described b y N emat-N asser and
Shok ooh ( 1979 ) led me t o in v es tig ate bo t h t he process t hat causes w ater pres-
sure chang es, and met hods of unders tanding w a v e dissipation. The process
of w a v e dissipation in soil is difficult t o model wit hout a consider able body
of ph ysical experiments designed t o inf or m t he model building process. Pro-
cesses t hat could contribute t o w a v e dissipation in sand include frictional
5
r ubbing of t he g r ains, and reflection, scattering, and g eometric spreading of
t he w a v e. Processes t hat could contribute t o reflection and scattering include
t he breakag e of f orce chains betw een particles t hat could contribute t o bo t h
t he production of sphericall y spreading high frequency har monics (due t o
im pacts betw een g r ains when f orce chains collapse) and increased frictional
heat g ener ation as g r ains r ub ag ains t each o t her . W it h such a wide v ariety of
processes t hat might contribute, a body of ph ysical experiments are required
in or der t o inf or m t he modeler .
Contr ibution T o Liq uef action The model de v eloped in chap ter 4 describes
ho w w ater press ure de v elops under an assumed g r ain def or mation, which
enters in ter ms of t he r ate of chang e of φ t he porosity o f t he material. This
chap ter giv es more detailed back g round on liquef action his t or y and t hen
proceeds t o de v elop an equation f or t he r ate of chang e of pore pressure
from firs t principles, full y incor por ating t he t her modynamics of w ater and
Darcy’ s la w , a w ell kno wn model f or viscousl y dominated flo w . In or der t o
close t he sys tem f or prediction from firs t principles, a model f or
@ φ
@ t
w ould
be needed, and t his model w ould require coupling a ph ysical descrip tion of
t he g r ain mo tion t o a continuum descrip tion of t he whole soil deposit.
Problems Wit h Continuum Models of Sand Ho w e v er , e v en if w e had a
good unders tanding of t he inter -g r ain inter actions, it is har d t o unders tand
ho w t o build a continuum model of g r ains inter acting in a w a y t hat mak es
consis tent ph ysical sense. W e w ant a continuum model, because a model
built on t he mechanics of individual g r ains depends on t oo man y details t o
be w ort hwhile. In particular , t o be accur ate w e need t o kno w t he irregular
shape, and material properties of an enor mous number of g r ains t hat fill
t housands of cubic meters of soil. A continuum model has t he property t hat
f or an y location in space, t here is some v alue of a rele v ant ph ysical v ariable,
such as density , v elocity , and so f ort h. In a netw or k of g r ains, it’ s clear l y
t he case t hat at a giv en point eit her t here is a g r ain, or t here is no g r ain,
so t hat if w e are t o proper l y describe t he sys tem wit h a continuum model
w e need t o unders tand ho w t he continuum model is t o correspond t o t his
discrete reality . In later chap ters of t his dissertation I de v elop a fr ame w or k
f or inter preting continuum models t hat can inf or m t he process of building
continuum models of discrete particles in a consis tent w a y .
6
Building an Exam ple Continuum Model of Discrete P ar ticles A model of
dissipation in t he sim ple Lennar d-Jones model of molecular inter actions is
an easier place t o s tart de v eloping t hese modeling techniques, and t he model
de v eloped in chap ter 5 of t his dissertation inf or ms us of an im portant ph ys-
ical process in w a v e dissipation, namel y momentum diffusion b y coupling
of t her mal and w a v e mo tion. This model w as built b y explicitl y using t he
techniques I de v eloped f or inter preting continuum models.
As a firs t pass at unders tanding t he effect of momentum diffusion, w e
can consider ho w kinetic ener gy t hat is tr ansf erred from a small mass t o a
lar g er mass ine vitabl y leads t o dissipation. Consider a 3g bullet tr a v eling at
500 m/s im pacting a 1 k g block of cla y on a frictionless table, and embed-
ding itself. The bulk -scale (non-t her mal) kinetic ener gy of t he sys tem bef ore
im pact is 0: 003 500
2
= 2 = 37 5J and t he momentum of t he sys tem is 0: 003
500 = 1: 5kgm=s . The momentum is conser v ed during im pact, so t hat t he
final v elocity of t he cla y wit h bullet embedded is ( 1: 5kgm=s)=( 1: 003kg) =
1: 496m=s and t he kinetic ener gy is 1: 003 1: 496
2
= 2 = 1: 12J or a loss of
( 1 1: 12= 37 5) 100 = 99: 7% . If momentum is conser v ed, but becomes
dis tributed o v er a lar g er mass, ine vitabl y t he kinetic ener gy associated wit h
center of mass mo tion mus t decrease, and t his is t he effect t hat predicts t he
dissipation in t he dissipativ e w a v e model, t hough t he causes are much dif-
f erent from t hose oper ating during a bullet im pact.
In t he next chap ter w e will visit t his model building technique as back -
g round f or t he models described in later chap ters.
7
Chapter 3
T echniques and Philosophical
Considerations f or Continuum M odels
If y ou do no t look at t hings on a lar g e scale it will be difficult f or y ou t o
mas t er s tr at egy
— Miy amo t o Musashi
It isn ’ t g ener all y t he pr actice of Engineers t o del v e t oo deepl y int o t he philo-
sophical im plications of t he models or procedures t he y use. Mos t of t he time,
t hese sorts of ques tions are sim pl y unnecessar y . W e can s tart wit h some ba-
sic descrip tion of phenomena which is kno wn t o w or k reasonabl y w ell in
pr actice, and w e can mak e modifications or predictions based on t hose de-
scrip tions t hat are inter preted im plicitl y in some concep tual fr ame w or k t hat
is common t o our colleagues. F or exam ple, t he N a vier -S t ok es equations f or
t he mo tion of a N e wt onian fluid are enor mousl y successful. W e can mos tl y
design air planes wit hout w orr ying about t he f act t hat air is actuall y a collec-
tion of individual molecules as t his f act has no measur able consequences at
t he scale of a ten meter long wing. Occasionall y t hough, w e will ask ques-
tions where a continuum descrip tion has ques tionable meaning. What is t he
meaning of a continuum descrip tion of t he flo w of a fluid t hrough a nano-
tube onl y 10 at omic diameters wide? Ho w does a continuum model of flu-
ids mak e sense at t he edg e of t he atmosphere where t he mean free pat h of
a molecule is per haps 1 meter? And especiall y in g r anular materials, ho w
do w e describe t he flo w of fluid and g r ains as a continuum when w e wish t o
describe f eatures on a similar lengt h scale as t he g r ains t hemsel v es?
8
3.1. Continuum Model as S tatis tical Mechanics
Jacob Bear explicitl y deals wit h t his ques tion in t he firs t chap ter of his
book on fluid flo w in porous media (Bear 1988 , p.20). T o aid in t he de v el-
opment of models in later chap ters, I will introduce here se v er al unifying
no tions t hat allo w us t o inter pret continuum models e v en when t he non-
continuum nature of t he material is apparent. In addition t o allo wing f or
better com parisons betw een reality and t he model predictions, t hese philo-
sophical and mat hematical issues can inf or m t he de v elopment of ne w mod-
els in t he initial s tag es.
3.1 Continuum Model as S t atistical Mechanics
Whene v er w e discuss a continuum model of ph ysical phenomena, w e neces-
saril y eit her ignore or mus t confront t he ques tion of what does a continuum
mean? Clear l y our moder n kno w ledg e tells us t hat air is made of near l y non-
inter acting molecules bouncing around wit hin t he entire v olume containing
t hem. W ater is made of molecules which inter act t o maintain a close pro x-
imity t o each o t her but do no t ha v e long-r ang e or der o v er man y molecular
diameters. The v olume t he y occup y s ta ys relativ el y cons tant as t he container
size increases due t o t he f or mation of a g as-liquid interf ace. A glass is lik e
a liquid in t hat it has molecules in close pro ximity and no long r ang e or der ,
but t he molecules are held in place in some kind of local or der which is pre-
ser v ed o v er long time scales. Ov er extremel y long time scales, in t he cr us t of
t he eart h f or exam ple, w e ma y ha v e material flo w in a material which w ould
nor mall y be considered a solid such as portions of t he eart h ’ s mantle. A
pol y cr ys talline solid has molecu les arr ang ed in a lattice wit h a cons tant lat-
tice orientation o v er a significant multiple of t he lattice spacing, but v ar ying
lattice arr ang ements from place t o place wit hin t he solid. A single cr ys tal
such as a quartz miner al ma y ha v e a cons tant lattice alignment t hroughout
its entire extent. These f or ms of matter mak e up t he bulk of e v er yda y mate-
rials, and y et w e frequentl y describe t heir beha vior in ter ms of a continuum.
The mat hematical property of a continuum is t hat it has a specific v alue f or
some v ariables of interes t at an y real number spatial coor dinate. Since in a
N e wt onian descrip tion, a molecular material can at bes t be said t o ha v e a
particular property onl y at exactl y t he coor dinates of t he molecules, clear l y
t he continuum is no t a representation of t he molecules t hemsel v es.
An alter nativ e approach is t o consider t he quantities of interes t in a con-
tinuum model as s tatis tical a v er ag es o v er some v olume t hat is small wit h
9
3.1. Continuum Model as S tatis tical Mechanics
respect t o t he o v er all body in ques tion but lar g e enough wit h respect t o t he
elementar y particles t hat mak e up t he body t o ha v e w ell controlled s tatis tics.
A continuum model of a ph ysical body is t hen a descrip tion of a lar g e but
unspecified number of t hese small v olumes. A v er ag es are t he appropriate
type of s tatis tic due t o t he presence of conser v ed quantities. The t o tal mass,
ener gy , momentum, and angular momentum of a sys tem under or dinar y
ph ysical conditions is conser v ed in t he absence of exter nal inter actions, and
t he linear nature of N e wt on ’ s La w s
∑
F = m a mak es t he t o tal of a conser v ed
quantity or its flux t he appropriate quantities of interes t. The property t hat
mak es a v er aging meaningful f or t hese types of quantities is:
N
∑
i= 1
x
i
=
K
∑
j= 1
n
j
⟨ x⟩
j
In o t her w or ds, t he t o tal o v er all N elementar y quantities x
i
of interes t,
when t hat t o tal is a conser v ed quantity , can also be calculated as t he sum
of all t he K extr apolated a v er ag es of g roups of v ar ying size n
j
. W e use t he
a v er ag e when w e are interes ted in t he beha vior of t he t o tal, so t hat s tatis tics
lik e t he median or o t her quantiles of t he dis tributions are no t appropriate in
t he presence of conser v ed quantities, ex cep t in so f ar as t he y ma y be effectiv e
appro ximations of t he a v er ag e.
This s tatis tical outlook on continuum models is no t ne w . The no tion of a
“representativ e v olume element” (R VE) is common in t he liter ature of multi-
scale modeling (cf. (Gal v anett o and Aliabadi 2010 )). Ho w e v er , it should be
no ted t hat t here is a continuum of sizes of v olume element, each wit h its
o wn deg ree of “representativ eness”. A tr ade off e xis ts betw een t he size of
s tatis tical fluctuations expected in t he s tatis tic, and t he size of t he smalles t
spatial detail t hat t he model can cap ture. This tr ade off can be illus tr ated
easil y in figure 3.1 . In t hese g r aphs, 10000 points w ere placed unif or ml y
at r andom wit hin t he inter v al ( 0; 1) . Because t he points are non-inter acting
t heir density is ex ceedingl y noisy com pared t o typical inter acting particles
which will tend t o arr ang e t hemsel v es int o a more or dered and less ener -
g etic s tate. The number density of points can be com puted b y taking a bo x
centered at a giv en location wit h v ar ying size. In addition t o t he unif or ml y
w eighted a v er ag e, it is also possible t o use a smoo t h k er nel w eighted a v er ag e
t o g et a lo w er noise s tatis tic, but f or illus tr ation pur poses t he unif or m bo x is
con v enient.
10
3.1. Continuum Model as S tatis tical Mechanics
9200
9400
9600
9800
10000
0.00 0.25 0.50 0.75 1.00
Box Size
n/size
Point Density At x = 0.5
9500
9750
10000
10250
10500
0.00 0.25 0.50 0.75 1.00
x
rho
Box Size
0.25
0.5
Density at a given point
F igure 3.1: As bo x size increases, t he fluctuations in density as a function
of bo x size decreases. In t he vicinity of size 0.25 t he density at x = 0.5 is
relativ el y insensitiv e t o t he specifics of t he bo x size (lef t g r aph). A t bo x size
0.25 t he density has a lar g e neg ativ e g r adient in t he vicinity of location 0.5
(right g r aph), f or bo x size 0.5 t his g r adient is smaller and relativ el y unif or m
across t he domain.
In an y mat hematical modeling procedure it is necessar y t o ignore certain
effects which are deemed t o be negligible, because if w e attem p t t o include
all details, our models become intr actable f or sys tems be y ond a f e w t hou-
sand or million molecules. Ev en in MD simulations po tentials are typicall y
tr uncated so t hat w e ignore v er y small long-r ang e f orces. The negligibility
of an effect is alw a ys wit h ref erence t o t he size of some o t her effect consid-
ered more im portant. In t he g r aph on t he lef t (figure 3.1 ), t he fluctuations in
t he density of t he material at x = 0: 5 when s 0: 25 is appro ximatel y 1= 10
of t he o v er all a v er ag e of 10000 points per unit of x , and t he a v er ag e v aries
on t he or der of 1= 100 as s chang es b y 10% or so. By t he time t he bo x size
is 0: 5 a v ariation of 10% in bo x size produces a v ariation in density of on
t he or der of 2% or so. F urt her more, t he v ariation at s = 0: 25 is related t o
t he f act t hat t here are more points slightl y t o t he lef t of 0: 5 t han t here are t o
t he right, producing a g r adient across t he domain. This sho w s t hat e v en t he
concep t of a “s tatis tical fluctuation” is tr aded off wit h t he concep t of spatial
inhomog eneity in t his process.
The continuum model t hen can be t hought of as an inter mediate asym p-
t o tic model (Barenblatt 2003 ) in which t he size of t he representativ e v olume
11
3.2. N on-dimensionalization, Multiple Scales, IS T and N ons tandar d
N umbers
is lar g e enough wit h respect t o t he minimum f eature size t hat s tatis tical
fluctuations in quantities of interes t are o( 1) relativ e t o t he typical v alues
of t he quantities, but t he R VE is small enough wit h respect t o t he o v er all
phenomenon of interes t t hat predicted results v ar y from place t o place or
time t o time at similar scales t o t he measurements which might be used t o
f alsify t he model. If no such inter mediate r ang e of scales exis ts, t hen a con-
tinuum model can no t be applied, and an explanation in ter ms of t he more
fundamental elemental particle descrip tion mus t be used.
It w ould be useful t o ha v e a mat hematical fr ame w or k in which w e could
express t he f act t hat t here is no exact size of t he “ideal” R VE. W e will fre-
quentl y w ant t o sim pl y assert t hat our “R VE” is “small enough” or t hat t he
error betw een measurements and t he predictions of our model are “small
enough” when t he number of R VE elements is “big enough. ” Such a con-
cep t mus t ultimatel y depend simultaneousl y on t he o v er all phenomenon of
interes t, t he measurement appar atus ’ size and precision, and t he deg ree of
s tatis tical uncertainty required b y t he user t o reject t he predictions of t he
model. This means t hat t here is no precise number t hat can be assigned uni-
v ersall y t o concep ts such as “small enough” or “big enough. ” In t he next
section I will introduce t he nons tandar d number sys tem de vised b y Edw ar d
N elson in his sys tem f or nons tandar d anal ysis kno wn as IS T (N elson 1977 )
and sho w ho w it can inf or m our modeling and v alidation procedures.
3.2 N on-dimensionalization, Multiple Scales, IST
and N onst andard N umber s
As pointed out b y man y aut hors, e v er y ph ysicall y meaningful equation mus t
be written as a sum of one or more ter ms in a homog eneous set of units
(cf. (Barenblatt 2003 ; F o w ler 1998 ; Ho wison 2005 ; Maha jan 2010 )). This
means t hat e v er y ph ysicall y meaningful equation can be re written in a di-
mensionless f or m b y dividing a dimensional equation b y a cons tant scale
f act or wit h t he same dimensions as t he ter ms. In particular , based on prior
kno w ledg e of t he relativ e sizes of t he ter ms in t he phenomenon of interes t
w e can usuall y find some meaningful combination of t he ph ysical v ariables
S , which has t he appropriate dimensions and f or which in typical condi-
tions max
i
j x
i
= Sj = O( 1) . That is t o sa y t hat an y ter m x
i
in t he equation is
under usual conditions about t he size of S or smaller . Af ter con v erting our
12
3.2. N on-dimensionalization, Multiple Scales, IS T and N ons tandar d
N umbers
equation t o non-dimensional f or m wit h appropriate scaling, w e can of ten
mak e prog ress b y asserting t hat certain ter ms ha v e an effect t hat is negligi-
ble and remo ving t hese ter ms from t he equation. N egligible, lik e beautiful,
is no t a w ell defined f or mal concep t, it mus t be wit h respect t o some willing-
ness t o mak e errors of some sort. Ho w e v er , w e can f or malize t he concep t of
negligible b y ref erence t o t he concep t of “infinitesimal” as a fr action of t he
ph ysicall y meaningful scale S .
Edw ar d N elson introduced a sim plified sys tem f or defining infinite (or
“unlimited”) and infinitesimal numbers in his paper (N elson 1977 ). This sys-
tem pro vides mos t of t he po w er of t he ear lier sys tem defined b y R obinson
( 1996 ) wit hout much of t he dizzying f or mal logic f oreign t o mos t applied
mat hematicians. That approach is elabor ated upon in application t o basic
mat hematical concep ts of calculus b y R obert ( 2011 ), and in modeling appli-
cations b y Lobr y and Sari ( 2008 ). This approach t o calculus melds w ell wit h
t he deriv ation process f or ph ysical equations especiall y ones in v ol ving mul-
tiple scales.
Useful Mat hematical Concepts from IST
Im portant concep ts from IS T are t he predicate “s tandar d” written st() . All
t he usual objects of nor mal mat hematics are “s tandar d” and t hose concep ts
defined wit hout ref erence t o t he predicate st() are called “inter nal”. An y
concep t which is defined explicitl y or im plicitl y using t he predicate st( x) is
called “exter nal”. Due t o t he w a y in which t hese concep ts are defined, w e
ha v e t he exis tence of a nons tandar d integ er , and can sho w t hat an y nons tan-
dar d integ er is bigg er t han an y s tandar d integ er , and hence is called “infi-
nite” or “unlimited” (N elson 1977 ). F urt her more, it f ollo w s t hat if N is a non-
s tandar d integ er t hen 1= N is a nons tandar d r ational closer t o zero t han an y
s tandar d number . It is numbers lik e t hese which are called “infinitesimal. ”
N umbers are t heref ore classified b y whet her t he y are “limited” or “unlim-
ited” and “infinitesimal” or “appreciable”. A number which differs from a
s tandar d number b y an infinitesimal amount is called “near -s tandar d”. Ev -
er y near s tandar d number x is infinitesimall y close t o a unique s tandar d real
number called t he “s tandar d part” of x and deno ted st( x) . In particular t he
s tandar d part of an infinitesimal number is 0 . These concep ts are defined
b y N elson ( 1977 ), and elabor ated in t he pre vious ref erences (Lobr y and Sari
2008 ; R obert 2011 ).
13
3.2. N on-dimensionalization, Multiple Scales, IS T and N ons tandar d
N umbers
The point o f such a number sys tem from our perspectiv e is t o tur n t he
usual calculus and asym p t o tic anal ysis int o an alg ebr aic manipulation of
special numbers and t hereb y f acilitate clear t hinking in de v eloping and an-
al yzing models. Since IS T has been pro v en t o be a conser v ativ e extension of
Zer melo-F r aenk el set t heor y wit h t he axiom of Choice (ZFC) t here is no t hing
ne w t hat can be pro v ed wit hin ZFC b y use of IS T (N elson 1977 ). Ho w e v er ,
t here ar e ne w met hods of proof which can f acilitate t he de v elopment of mod-
els. Jus t as t he C prog r amming languag e is com piled t o machine code and
t heref ore is no t able t o represent an y ne w com puter prog r am t han t he ones
representable b y writing r a w machine code, w e ne v ert heless find t hat C is
a more con v enient languag e in which t o prog r am t han writing r a w binar y
machine code because it mak es concep ts clearer t o t he humans r eading t he code .
The “effectiv e infinitesimal” and modeling er ror s
An “infinitesimal” number in IS T is f or mall y so small as t o be indis tinguish-
able from zero wit hout t he fine-g r ained comb of t he st predicate so t hat no
mat hematician w or king in a s tandar d mat hematical fr ame w or k could dis tin-
guish it from zero. Ho w e v er , in pr actical applied ter ms, it is a f or mal logical
pro xy f or t he concep t of negligible wit hin t he scale of interes t defined b y t he
modeler . Declaring sim pl y t hat w e will treat an expression as infinitesimal
allo w s us t o t hen deter mine alg ebr aicall y what o t he r expressions could also
be considered as infinitesimal. It defines in some sense a separ ation of scales
int o t hose w e can obser v e or care about, and t hose t oo small t o matter f or our
pur poses. When taking t his approach, w e mus t t heref ore be careful when
multipl ying quantities ε w e ha v e arbitr aril y tagg ed as “effectiv el y infinitesi-
mal” b y v er y lar g e numbers of or der O( 1= ε) since t hese mus t be considered
as “effectiv el y infinite”.
Since w e will alw a ys tr ansition from t he languag e of IS T ultimatel y back
t o t he languag e of s tandar d mat hematics, t he process of mo ving int o IS T
(called tr ansf er ), de v elopment of a model, and mo ving back t o s tandar d mat h-
ematics (called s t andar dization ) will if w e are talking about ph ysical models
where no quantities are actuall y infinitesimal or infinite, aut omaticall y mean
t hat w e are committing errors. The goal of an y mat hematical modeling ef-
f ort is t o bound t he extent of t hese errors t o be small enough t hat treating
t hem as if t he y w ere infinitesimal is pr acticall y jus tified.
One no tion t hat t he infinite numbers of IS T can help us wit h is t he near -
com plete separ ation of different ph ysical scales as represented b y t he relativ e
14
3.2. N on-dimensionalization, Multiple Scales, IS T and N ons tandar d
N umbers
size of nons tandar d integ ers. If N is a nons tandar d integ er , and M is a differ -
ent nons tandar d integ er , t hen N= M might be infinitesimal, appreciable, or
infinite depending on t he relativ e sizes of t hese quantities. If M = N
2
t hen
N= M = 1= N 0 whereas if M =
p
N t hen N= M
p
N 1 , and finall y
if M = aN where a is a s tandar d number , t hen N= M = 1= a a perf ectl y good
appreciable s tandar d number .
Exam ples of Scale Separation
F or exam ple, w e ma y be interes ted in modeling t he mixing of t he ocean. The
mean dep t h of t he ocean ma y become our lengt h scale of interes t so t hat one
non-dimensional unit of lengt h represents about 4km. If w e wish t o repre-
sent t he dissol ving of carbon dio xide at t he atmospheric interf ace, w e ma y
be interes ted in a lengt h scale of onl y a f e w tens of at oms 1 nm. This
is appro ximatel y 2: 5 10
13
non-dimensional units. T reating t his lengt h as
infinitesimal will in v ol v e essentiall y no meaningful error in our model. If
w e wish t o include t he s tirring and mixing effect of s trong wind and w a v es
w e ma y need t o a v er ag e o v er per haps up t o 10 m w a v e cres t height. This is
0: 0025 non-dimensional units, a quantity which w e ma y choose t o also treat
as infinitesimal on t he o v er all dep t h scale. F inall y , per haps bubbles f or med
during t he breaking of w a v es are typicall y 1cm in diameter but can be mix ed
t o a dep t h of t he 1 0m w a v e height. The processes occurring in t he vicinity
of a single bubble tak es place on a scale of per haps 10cm or 2: 5 10
5
non-
dimensional lengt h units and a single bubble might be considered infinites-
imal relativ e t o t he 10m w a v e height. In t his situation, t hese infinitesimal
scales, t he molecular mixing, bubble inter action, and w a v e mixing scale are
extremel y dispar ate so t hat each ma y be considered infinitesimal relativ e t o
t he next lar g es t. Hence w e ha v e dep t hs of or der 1 betw een t he surf ace and
t he bo tt om of t he ocean, of or der ε 0 near t h e surf ace where w a v es break , of
or der appro ximatel y ε
2
0 at t he bubble scale, and of or der appro ximatel y
ε
5
0 at t he molecular scale.
Ano t her im portant concep t f or modeling t hat w e g ain from IS T is t he
concep t of “s-continuous” (roughl y meaning “seems continuous”). An s-
continuous function chang es b y onl y an infinitesimal amount when its in put
is chang ed b y an infinitesimal amount (R obert 2011 ). Ho w e v er , especiall y
f or our pur poses, t his ma y be tr ue of some non-s tandar d functions t hat are
no t actuall y continuous, such as a s tep function which is cons tant betw een
infinitesimal s teps, but t he size and spacing of t he s teps is at an infinitesimal
15
3.2. N on-dimensionalization, Multiple Scales, IS T and N ons tandar d
N umbers
0.5
1.0
1.5
2.0
2.5
0.00 0.25 0.50 0.75 1.00
X
Y
A continuous function and
an s−continuous type approximation
F igure 3.2: An appro ximation of a continuous function (black) b y fine
g r ained s teps (red).
scale. P er haps t he difference betw een t he s-continuous function and some
s tandar d continuous function remains infinitesimal at e v er y point. W e can
appro ximate t his concep t in t he g r aph of figure 3.2 . In f act, t he difference
betw een t he s tep function and t he continuous function has been accentuated
here t o mak e it more visible com pared t o an ear lier v ersion where t he s tep
sizes w ere smaller b y a f act or of 5.
Mat hematicall y , each bounded s-continuous s tep-wise non-s tandar d ap-
pro ximation corresponds t o some unique s tandar d function t hat is continu-
ous since e v er y s-continuous s t andar d function is continuous (R obert 2011 ).
In t he context of t he deriv ation of ph ysical continuum models, a small bo x
in t he vicinity of a point defines t he region o v er which an a v er ag ed ph ysical
quantity is defined. F or t he pur poses of modeling w e ma y treat t he quan-
tity defined wit hin t he bo x as a cons tant o v er t he bo x. The difference be-
tw een adjacent bo x es when divided b y t he displacement betw een t he cen-
ters of t he bo x es can define a deriv ativ e up t o an infinitesimal error . The
difference betw een t his deriv ativ e of an s-continuous s tep function and t he
deriv ativ e of t he actuall y continuous function t hat is t he s tandar dization of
t he s-continuous function can onl y be dis tinguished mat hematicall y wit hin
16
3.3. The Model Cons tr uction T em plate: Measurement and Molecules
t he context of IS T .
In building continuum models as models f or s tatis tical a v er ag es, w e will
consider small bo x es of nor mall y unspecified size defining t he spatial ex-
tent of t he a v er ag ed quantities, deriv e equations f or t he e v olution of t he a v -
er ag e quantities using t he discontinuous s tep-function appro ximation im-
plied, treating t hese functions as s-continuous, and t hen con v ert t he result-
ing equations int o t heir corresponding s tandar d equations b y t he process of
s tandar dization.
Since real ph ysical quantities are ne v er actuall y infinitesimal, t he entire
process will im pl y spatial, time, or ener gy scales belo w which t he equations
are no t expected t o mak e exact predictions. If t he sizes of our errors are
lar g e enough t hat w e wish t o account f or t hem explicitl y , a s t ochas tic com-
ponent can be included which is infinitesimal onl y on a v er ag e but whose
error at e v er y point can be described b y r andom v ariables t hat are appre-
ciable and po tentiall y correlated. Thus, when describing a ph ysical sys tem
using a continuum model t here is ine vitabl y a sub-scale s tr ucture “wit hin
t he bo x” about which w e onl y ha v e a v er ag e inf or mation. There are man y
po tential sub-scale s tates t hat are consis tent wit h t he a v er ag e and are natu-
r all y described as if t he y w ere r andom since w e ha v e little t o no inf or mation
wit h which t o predict t hem. Depending on t he specifics of t he sys tem, a v ari-
ety of dis tributions and correlation s tr uctures ma y be applicable t o describe
t hese fluctuations from t he predicted v alues.
3.3 The Model Constr uction T em plate:
Measurement and Molecules
The Measurement R egion Size
W e begin our model cons tr uction tem plate wit h t he no tion of a control v ol-
ume (CV) or representativ e v olume element (R VE). Ho w e v er , w e specificall y
inter pret t his R VE as being t he same or der of magnitude in size as t he region
o v er which our t heoretical or actual measurement appar atus a v er ag es. In
an y giv en scientific problem, w e can onl y v alidate or f alsify our model b y ref-
erence t o ph ysical measurements, a f act sometimes ignored in purel y mat h-
ematical anal ysis of continuum models. But t he appropriate ph ysical mea-
surement appar atus is necessaril y different from one problem t o ano t her .
A non-contact t her mometer f or exam ple measures t he t o tal infr ared r adia-
17
3.3. The Model Cons tr uction T em plate: Measurement and Molecules
tion com ing from some region of a plate, a pressure meter measures t he t o tal
f orce on some small flexible membr ane of kno wn area, a digital camer a tak es
imag es whose resolution is limited b y t he pix el size and t he op tical proper -
ties of t he lens, in each case, a measurement is an integ r ation procedure,
a v er aging o v er some region of space and time. In t he fines t possible mea-
surements w e ma y g et measurements of t he appro ximate locations or f orces
acting on individual at oms or molecules so t hat our model mus t produce
predictions at t he individual at omic le v el . Ho w e v er , in t he common case,
t he spatial and time scale o v er which measurements occur is much lar g er
t han t hose rele v ant t o an at om. Alt hough w e of ten describe continuum mod-
els as “infinite dimensional”, t he real v alue of a continuum model is t hat it is
a tem plate f or or f amil y of finite dimensional models which can giv e accur ate
predictions e v en f or t hose members of t he f amil y whose number of dimen-
sions is v as tl y smaller t han t he dimensionality of t he detailed MD simulation
being appro ximated. F urt her more t he con v er g ence of t he finite dimensional
models as dimension increases im plies t hat answ ers are appro ximatel y inde-
pendent of t he number of dimensions f or some lar g e enough dimensionality .
Dimensionless Lengt h R atios
Once w e define a measurement region size, t here are no w at leas t tw o lengt h
scale r atios of interes t in t he problem. The firs t is l
m
= l
b
t he size of t he mea-
surement region com pared t o t he o v er all body of interes t.
The next is l
a
= l
m
t he r atio of t he inter at omic or inter -particle dis tances t o
t he size of t he measurement region. This r atio also deter mines t he typical
number of elements in an R VE N/ l
3
m
= l
3
a
f or a 3D problem or l
2
m
= l
2
a
f or a 2D
problem. It should be no ted t hat w e ma y ins tead t hink of l
a
as t he “inter -
particle” dis tance if w e are considering a model f or discrete particles lik e
pellets or sand g r ains or inter acting colloidal suspensions.
As an exam ple of some of t hese r atios in pr actice, suppose w e are in-
teres ted in eart hquak e w a v es. W e ha v e measurements via a seismometer at
100Hz sam pling frequency . The Nyquis t frequency is t heref ore 50Hz , and at
around 3000m=s w a v e speed, one v elocity measurement represents t he a v -
er ag e v elocity o v er a spatial region of or der 60m . Clear l y 60m is enor mous
com pared t o t he inter at omic dis tance in rock , y et 60m is tin y com pared t o
t he o v er all lengt h of t he w a v e tr ain las ting per haps 30s and t hereb y extend-
ing o v er 90000m or t o per haps 10
6
m t hat t he w a v e might tr a v el t hrough t he
eart h ’ s cr us t t o our seismometer . On t he o t her hand, f or a small mechanical
18
3.3. The Model Cons tr uction T em plate: Measurement and Molecules
Ratio of Measurement Size to Atomic Element Size
Ratio of Body Size to Measurement Size
F igure 3.3: V arious regimes of size r atio deter mine which types of models
are appropriate. In t he upper diag r ams t he bo x represents t he measurement
size, and t he circles represent molecules or discrete elements. In t he lo w er
diag r ams t he bo x represents t he body and t he g re y square represents t he
region o v er which measurements a v er ag e.
part, per haps a s tr ain g aug e is on t he or der of 1mm in size f or a part on t he
or der of 100mm so t hat t his r atio is more moder ate. Sometimes t he mea-
surement size is essentiall y t he same size as t he sys tem. F or exam ple, a lar g e
pol ymer biomolecule wit h a super -fine laser measurement, or a pis t on in a
machine wit h a single f orce g aug e on t he rod.
F inall y , t here is t he ques tion of t he inter action lengt h scale l
i
= l
a
. T ypicall y
t his r atio is > 1 but can sometimes be v er y lar g e. When t his r atio is small,
onl y a f e w neighbors of an at om contribute t o t he f orces on a molecule, but
when it is lar g e, such as when significant electromagnetic f orces are in v ol v ed
on bare ions, or when w e model t he g r a vitational inter actions of as teroid
clouds or s tar clus ters, t he t o tal f orce on a particle mus t include contributions
from f ar a w a y from t he particle.
In addition t o lengt h scale r atios, t here are also timescale r atios of in-
teres t. An im portant one is dt
meas
= t
where dt
meas
is t he time o v er which a
measurement de vices a v er ag es when it giv es us a measurement, and t
is
t he t o tal time t hat w e obser v e t he sys tem and wish t o mak e predictions or
explain t he sys tem wit h t he model. Moder n measurement equipment w or k s
19
3.3. The Model Cons tr uction T em plate: Measurement and Molecules
b y filtering analog electrical signals t o bandwidt h limit t hem, and t hen sam-
pling t his bandwidt h limited signal. The filtering inherentl y delocalizes t he
signal in time b y an amount on t he or der of 1= f
cutoff
and sam pling induces
some additional issues including jitter and con v ersion time, so t hat e v en in a
classical setting w e are subject t o Heisenber g type uncertainty betw een fre-
quency com ponents and t heir location in t he tim eseries as is w ell kno wn.
Analog measurement techniques do no t pro vide an y w a y around t his issue
as w e mus t s till someho w define what w e mean b y t he exact v alue of t he
measurement and at what time t he measurement occurred.
In t he subsequent sections, whene v er a lengt h is mentioned it is tak en
t o mean a r atio of t he related dimensional lengt h t o some o v er all lengt h of
t he lar g e scale sys tem l
b
, and whene v er a time is me ntioned it is tak en t o
mean as a fr action of t he whole time of some experiment t
. These r atios
are t heref ore dimensionless quantities and g ener all y all less t han 1 in size.
Occasionall y w e will mak e t his explicit in s tatements such as l
m
= l
b
≪ 1 which
can be triviall y inter preted as ( l
m
= l
b
)=( l
b
= l
b
) ≪ 1 where s tarred quantities
are dimensional quantities in some units of measurement such as SI.
The molecular nature of t he tr ue dynamics, and t he
appro ximate spatial continuity of t he measurement
N ext, w e consider t he N e wt onian equations of mo tion f or t he molecules in
some ph ysical sys tem of interes t (which w e assume ha v e cons tant mass):
dp
dt
= m
dv
dt
=
N
∑
i= 1
F
i
+ F
b
(3.1)
Here N ref ers t o all t he molecules wit hin a dis tance of t he or der of l
i
,
t he inter action l engt h and F
b
is a f orce t hat comes from regions f ar lar g er
t han t he inter action lengt h and is cons tant o v er t he measurement v olume,
g ener all y t his in v ol v es g r a vitational f orces or occasionall y electromagnetic
w a v es propag ating from f ar a w a y .
F inall y , w e consider t he “equation of mo tion” f or t he measurement of t o-
tal momentum P o v er K molecules in t he region of size l
m
. The linearity of
N e wt on ’ s la w s, and t he conser v ation la w s of momentum and ener gy con-
vince us t hat t he measurement mus t be inter preted as a v er ag es, since t he
a v er ag e can be linear l y extr apolated t o t he t o tal.
20
3.3. The Model Cons tr uction T em plate: Measurement and Molecules
dP
dt
=
d
dt
( K( l
m
)⟨ m v⟩) =
d
dt
0
@
K( l
m
)
∑
i= 1
m
i
v
i
1
A
=
K( l
m
)
∑
i= 1
0
@
N
i
( l
i
)
∑
j= 1
F
i j
1
A
+ F
b
= F
net
(3.2)
Con v er sion of t he discrete sum o v er indices to a nonst andard
integ ral o v er space
So f ar , t his dynamics tak es place in an abs tr act discrete t opology . Each parti-
cle is its o wn “open set” and is disconnected from t he o t her particles. Thank s
t o finite a v ailable ener gy , and our assum p tion of appro ximating quantum
mechanics t hrough N e wt on ’ s mechanics and an appropriate choice of inter -
action po tential, and assuming no nuclear fusion can occur , tw o particles
ne v er approach arbitr aril y close t o each o t her . Hence, each particle can be
modeled as a v er y tin y ball or bo x in 3D space. W e can write a mass dis tri-
bution as
ρ
discrete
( x) =
N
∑
i= 1
m
i
B( x x
i
) (3.3)
Where B = 1= ε
3
when( x x
i
) is inside t he cube of side ε around 0 so t hat
t he t o tal integ r al is 1. N o w , in IS T , w e can assume t hat ε 0 and reco v er t he
Dir ac delta function as a perf ectl y nor mal but nons tandar d function defined
pointwise. Ph ysicall y , no t hing is wrong wit h t his because w e do no t ha v e
enough ener gy t o bring tw o particles close enough t og et her t o disco v er t he
specific size of ε so t hat a nons tandar d infinitesimal is a v er y v alid model
of t he obser v able ph ysics. The momentum dis tribution can be defined sim-
ilar l y , b y associating a v ect or v elocity t o each particle.
P
discrete
( x) =
N
∑
i= 1
m
i
v
i
B( x x
i
)
These nons tandar d mass dis tribution functions are triviall y integ r able
pro vided t hat w e use a small enough elemental v olume, one whose side is
infinitesimal relativ e t o ε . Since t he Dir ac function is a nons tandar d function,
it is no t sur prising t hat t he integ r al depends on t he particular nons tandar d
spatial s tep size.
21
3.3. The Model Cons tr uction T em plate: Measurement and Molecules
Ho w e v er , t he scientific quantity of interes t is no t usuall y t he unmeasur -
able individual particle dynamics, but r at her t he dynamics of t he measure-
ment
dp
dt
+
dp
dt
transport
=
d
dt
0
@
∫
D
p
discrete
( x) dV( x)
ε
1
A
=
∫
x2 D
0
B
@
∫
y= x+ O( l
i
)
F( x; y) dV( y)
ε
1
C
A
dV( x)
ε
+ F
b
( x)
=
K( l
m
)
∑
i= 1
0
@
N
i
( l
i
)
∑
j= 1
F
i j
1
A
+ F
b
= F
net
(3.4)
This equation is v alid so long as cdt≫ l
b
which is t o sa y t hat an increment
of time which is small com pared t o t he t o tal obser v ation time, light can tr a v el
a dis tance which is much lar g er t han t he o v er all size of t he body , so t hat t he
retar dation of f orces, and o t her relativis tic effects can be ignored.
So f ar t he equations are exact
dP
dt transport
represents t he chang e in t he mea-
surement purel y caused b y molecules entering or lea ving t he measurement
region, in situations where vdt ≪ l
m
, t his might be neglected as molecules
drif t onl y a tin y fr action of a measurement dis tance in our infinitesimal time.
W e ha v e sim pl y replaced a s tandar d discrete sum o v er particle index es (eq
3.1 ), wit h a finite but nons tandar d sum (an integ r al) o v er locations in space,
b y associating each particle and its mass, momentum and pair wise f orce
inter actions wit h nons tandar d functions of space. N o te t hat if w e require
multi-body po tentials t o accur atel y appro ximate quantum effects, w e can
create a f orce k er nel t hat is a function of an y number of locations. F or sim-
plicity of exposition w e assume pair wise f orces are sufficient.
Ho w e v er , in t he absence of enor mous com puting po w er , and extr aor di-
naril y fine measurements f or initial conditions, t he required molecular dy -
namics calculations necessar y t o carr y out t he exact solution are com pletel y
prohibitiv e. Consider t he eart hquak e dynamics problem, where w e mus t
com pute t he molecular dynamics of t he entire eart h! A com puter t o do t his
is necessaril y going t o ha v e appro ximatel y t he same mass as t he eart h. F or -
tunatel y , t hese detailed calculations are usuall y unecessar y as w ell since in
22
3.4. The s tatis tical nature of t he model, f alsifiability , par ameter es timation,
and goodness of fit
t he absence of enor mousl y detailed measurements almos t all of t he com pu-
tational output is non-f alsifiable.
3.4 The st atistical nature of t he model,
f alsifiability , parameter estimation, and
goodness of fit
It goes wit hout sa ying t hat t o model an actual ph ysical experiment w e mus t
ha v e a measurement of some initial conditions, and of boundar y or f orcing
conditions t hat occur t hroughout t he experiment. T o v alidate t he model pre-
diction w e mus t ha v e some measured data t o com pare wit h predictions at
v arious space and time points. In t he model cons tr uction tem plate so f ar w e
ha v e em phasized t hat t he predictions of a continuum model are predictions
of spatial a v er ag es o v er small regions of space. As is t he case in a v er ag-
ing t here is alw a ys some de viation betw een t he a v er ag e v alue of a quantity
and t he particular v alue at some point. F ailure t o predict t he exact measured
quantities at an y giv en time point can no t b y itself be cons tr ued as a f ailure of
t he model. Ins tead, t o v alidate t he model w e mus t com pare t he predictions
and measurements relativ e t o some measure of ho w lik el y t h e de viation is
under some reasonable probabilis tic model f or t he v ariations, and also more
g ener all y under some usefulness criterion f or predictions wit h de viations of
t hat magnitude (a utility or decision model).
F requentl y some quantities in t he model are no t precisel y kno wn a-priori,
quantities such as t he Y oung’ s modulus of a material, t he precise lengt h of a
s tr uctur al element, t he char acteris tic time of some inter nal process, or t he ac-
tiv ation ener gy f or some chemical reaction in t he presence of a catal ys t. Also,
initial and f orcing conditions ma y no t be precisel y kno wn, or onl y kno wn at
a small number of points. F requentl y w e are in t he position of needing t o
es timate t hese unkno wn quantities from t he data and t hen use t hese es ti-
mated quantities t o predict what w ould or did happen under some condi-
tions of interes t. Alt hough t here are se v er al approaches t o t he application
of probability t heor y t o s tatis tical inf erence, one which is especiall y useful
f or ph ysical modeling is t he Ba y esian approach in which uncertainty about
a quantity which ma y be eit her v ariable or fix ed and un v ar ying from exper -
iment t o experiment can ne v ert heless be giv en a probability associated t o
different v alues, and inter preted as t he relativ e reasonableness of t he v ari-
23
3.4. The s tatis tical nature of t he model, f alsifiability , par ameter es timation,
and goodness of fit
ous v alues t hat t his quantity might ha v e actuall y tak en on.
V alidation in idealized detail under t his fr ame w or k is t hen t he process
of collecting data, in putting t hat data int o t he model’ s initial, boundar y , and
f orcing conditions, inf erring t he dis tribution of im portant quantities of in-
teres t, and t hen com paring t he predictions of t he model under high proba-
bility v alues f or t he unkno wn quantities t o t he data collected at v arious time
points, and assessing whet her t he de viations betw een t he predictions and
t he data v alues indicate a w ell calibr ated model or if t here is some sys tem-
atic f ailure of t he model t o predict w ell in certain regimes. Such sys tematic
errors or biases w ould indicate t hat a fundamental ph ysical process ma y be
lef t out or poor l y described in t he model. F inall y w e mus t deter mine whet her
an y apparent f ailures are of pr actical interes t, or if t he model predicts w ell
enough f or t he pur poses f or which it will be used.
One f act of im portance is t hat measurements are ne v er absolutel y pre-
cise. The fines t quality analog t o digital con v erters t hese da ys use on t he
or der of 30 or so bits of precision, but accur acy of t he lo w or der bits is onl y
guar anteed if t her mal and electrical noise are carefull y controlled in ins tr u-
ment cons tr uction. In an y case, w e can no t imagine a situation in which an y
measurement will e v er in t he his t or y or future of humans ha v e an ywhere
near 190 bits of precision, which corresponds t o counting t he pro t ons in t he
sun t o wit hin plus or minus one pro t on or so. It’ s trivial t o sa y t hen t hat s ta-
tis tical error mus t logicall y enter int o e v er y calculation, e v en if in t he end it is
of no pr actical consequence and a decision is t hen made t o explicitl y ignore
it. F rom t he Ba y esian per pectiv e it is possible t o assign probabilities t o mod-
eling error as w ell, so t hat t he sum of measurement and modeling error ma y
sometimes be lum ped int o a single quantity t hat can inf or m t he lik elihood
of t he data P( Dj a) where D is t he measured v alues of t he data, and a are
t he kno wn or unkno wn par ameters in t he model including t he v alues of ini-
tial, boundar y , and f orcing conditions, as w ell as t he uncertain fundamental
quantities mentioned abo v e.
It is common in almos t all circums tances t o drop ter ms from an equation
when t hose ter ms can be a-priori deter mined t o be of much smaller or der of
magnitude t han t he main ter ms. F or exam ple w e ma y ignore t he v ariability
in t he g r a vitational acceler ation betw een t he bo tt om of a building and t he
t op of a building when deriving an equation of mo tion f or objects f alling off
of a skyscr aper . Such appro ximations necessaril y define bo t h a bias and an
im plicit “minimum scale” f or allo w able de v iations betw een t he model and
t he data e v en if an extremel y precise nanosecond accur ate clock is used t o
24
3.5. Conclusion
time t he f all so t hat t he error in t he measurement ins tr ument w ould seem t o
be of smaller scale t han t he obser v ed de viations. Ex cep t in cases rele v ant t o
metrology or spacecr af t orbits or o t her areas where precision is absolutel y
necessar y , it is r arel y of pr actical benefit t o model a process t o accur acy better
t han about 1 part in 10
4
giv en t he difficulty of collecting sufficientl y accur ate
in put data f or such pur poses. This t hought process sho w s t hat alt hough w e
idealize certain quantities as “infinitesimal” in our use of IS T , in f act e v en
t hose who do no t adop t t he NS A approach already treat quantities whose
or der of magnitude is a perf ectl y s tandar d number as if t he y w ere infinites-
imal, tr uncating t hem out of t heir equations bef ore use.
3.5 Conclusion
The goal of t his chap ter w as t o introduce an inter pretation of continuum
models as models f or local spatial a v er ag es of elementar y properties, t o ar -
gue explicitl y f or t he application of t he a v er ag e due t o t he presence of con-
ser v ed quantities, and t o introduce t he concep ts and ter minology of IS T t o
giv e us t ools f or cons tr ucting and describing models t hat better match mat h-
ematical concep ts wit h ph ysical concep ts. In f ollo wing chap ters t hese con-
cep ts will be used as needed t o describe no v el models f or t he flo w of w a-
ter during soil liquef action, and f or t he dissipation of w a v es in a simulated
molecular sys tem which conser v es ener gy t o extremel y high precision. The
molecular dissipation model is meant t o be sugg es tiv e of met hods t hat could
ultimatel y lead t o a continuum model of ho w w a v e ener gy induces porosity
chang es, which in tur n couples t o, and induces pore pressure chang es and
fluid flo w during liquef action.
25
Chapter 4
The D egree of Drainage During
Liquef action
I am t he daught er of Eart h and W at er ,
And t he nur sling of t he Sky ;
I pass t hr ough t he por es of t he ocean and shor es;
I c hang e, but I canno t die.
— P ercy Bysshe Shelle y
In t he textbook definitions of soil liquef action, one assum p tion is par amount:
t hat during t he eart hquak e no significant quantity of w ater can flo w (Kr amer
1996 ; Holtz and K o v acs 1981 ). Ho w e v er , w ater pressure is a function of den-
sity and tem per ature (NIS T 2012 ). In particular , a linear T a y lor series f or
small de viations from some particular typical density and tem per ature al-
lo w s us t o predict t he chang es in pressure associated wit h liquef action as a
sim ple linear function of density and tem per ature. Ear l y in m y research I
t ook t he undr ained h ypo t hesis as f act since it is repeated widel y and con-
sis tentl y wit hin t he liter ature. Ho w e v er , in or der f or pressure t o chang e,
in t he absence of w ater flo w , w e mus t ha v e eit her t he g r ains g etting bigg er
iso tropicall y so t hat t he v olume of v oids decreases unif or ml y e v er ywhere,
or w e mus t ha v e heating of t he w ater . The firs t h ypo t hesis is untenable as
t here is no ph ysical basis f or such an effect. Heating, on t he o t her hand, is
tenable and I spent a consider able time t hinking about ho w t his effect might
w or k. Ev entuall y , ho w e v er , I deter mined t hat t he bes t s tarting point w as t o
deriv e an equation f or t he r ate of chang e of w ater pressure and deter mine
26
4.1. Back g round
which v ariables w ere actuall y responsible, including all effects t hat I could
remo tel y consider reasonable.
The s tarting point f or t his deriv ation w ere t he concep t of conser v ation of
mass, Darcy’ s la w , and t he T a y lor series f or w ater’ s pressure v s density and
tem per ature function.
The resulting model ultimatel y led t o t he f ollo wing paper which is repro-
duced here in t he f or m in which it w as submitted f or publication in Proceed-
ings of t he R o y al Society A (Lak eland, R echenmacher , and Ghanem n.d. ).
4.1 Backg round
The g ener al phenomenon of t he liquef action of g r anular materials is applica-
ble t o a wide v ariety of g r ains and fluids. Here, w e f ocus on one of t he mos t
im portant aspects of t his phenomenon, namel y t he role of w ater flo w dur -
ing eart hquak e induced liquef action of sand near t he g round surf ace. Our
anal ysis uses a nondimensional f or mulation and asym p t o tic anal ysis which
is adap ted t o t he firs t f e w tens of meters of soil where t he inter pla y of g r a v -
ity , g r ain rearr ang ement, w ater com pressibility , and w ater flo w combine t o
cause t hese des tr uctiv e e v ents.
In t he Geo technical Engineering field, soil liquef action is commonl y un-
ders t ood as a consequence of w ater pressure buildup due t o r apid squeezing
of pore spaces, wit hout sufficient time f or w ater t o flo w t hrough t he g r ains
and dr ain t he pressure, e.g. (Sa wicki and Mierczyński 2006 ; Kr amer 1996 ;
Holtz and K o v acs 1981 ). When g r ains are loosel y pack ed, during eart hquak e
mo tion t he tendency is f or soil g r ains t o mo v e closer t og et her , squeezing t he
w ater and r apidl y increasing t he pressur e due t o t he high bulk modulus of
w ater .
T ak en across a t hin horizontal section of t he soil, t he t o tal v ertical f orce is
t he sum of t he contact f orces betw een g r ains, and t he w ater pressure times
t he cross sectional area. This “t o tal s tress” minus t he w ater pressure is t he so
called “effectiv e s tress”, commonl y used in g eo technical anal ysis (Holtz and
K o v acs 1981 ), which is a measure of t he contribution of g r ain-g r ain inter ac-
tions in t he soil. If t he contact f orces drop t o zero, t hen t he w ater pressure
carries t he entire nor mal s tress, and shear s tresses induce flo w in t he manner
of a viscous fluid.
The goal of liquef action assessment up t o no w has been t o deter mine
ho w t he w ater pressure in t he soil will chang e during cy clic loading. The
27
4.1. Back g round
met hodology em plo y ed has primaril y been t o build models based on t he re-
sults of labor at or y triaxial, hollo w cy lindrical, and sim ple shear tes ts, as w ell
as more extensiv e ph ysical centrifug e models. Labor at or y tes ts use sam ples
of sand typicall y around 10 t o 20 cm in char acteris tic size. The sand is sur -
rounded b y an im per meable, flexible membr ane t o tr ap t he pore w ater , and
t he entire sam ple is contained in a pressurized v essel t o simulate t he o v er -
all pressure conditions in t he g round. Cy clic loading of v arious f or ms is
applied and t he t o tal and effectiv e s tress s tates, and e v olution of pore w a-
ter pressure are tr ack ed t hrough t he cy cling. These experiments ha v e been
extensiv el y perf or med. T o com plement t hese small labor at or y scale tes ts,
extensiv e ph ysical modeling in v arious types of g eo technical centrifug e ap-
par atus ha v e been carried o ut. These tes ts use centrifug al acceler ation t o
model t he s tresses induced in deposits of soil t hat are 10 t o 100 times deeper
t han t he model scale and allo w full 3 dimensional g eometries t o be simulated
wit hout t he ambiguity of numerical simulations of com plex soil materials.
Combined wit h t hese ph ysical obser v ations, engineer s ha v e em plo y ed cor -
relations wit h obser v ed field conditions in pos t-seismic in v es tig ations. An
o v er vie w of t he current s tate of liquef action research ma y be f ound in Sa w -
icki and Mierczyński ( 2006 ).
Im plicit in t he use of small labor at or y sam ples wit h im per meable mem-
br anes, and explicit in man y textbook definitions of liquef action, is t he as-
sum p tion t hat during t he eart hquak e t here is no significant w ater flo w or
chang e in w ater v olume, which is kno wn as t he “undr ained” condition (cf.
(Kr amer 1996 ; Holtz and K o v acs 1981 )). The assum p tion t hat soil liquef ac-
tion occurs under undr ained conditions has lead t o extensiv e research int o
undr ained tablet op experiments such as triaxial and sim ple shear tes ts, wit h
se v er al met hods sugg es ted t o o v ercome t he small v olume chang es allo w ed
b y t he com pliance of t he r ubber membr ane (Siv at ha y alan and V aid 1998 ).
Alt hough centrifug e and labor at or y tes t data ha v e long sho wn t hat w ater
mig r ation can occur (F ieg el and K utter 1994 ; K okusho 1999 ), t hese situa-
tions ha v e been treated as if t he y w ere ex cep tions t o t he nor mal situation
of undr ained pore pressure increase. Ho w e v er , giv en t he t her modynamics
of w ater pressure, t he undr ained assum p tion can lead t o ph ysicall y incorrect
predictions. If t he w ater is no t allo w ed t o chang e v olume whatsoe v er , t hen
t he onl y mechanism f or pressurization is heating. This sho w s t hat some care
is required t o deter mine proper l y t he role of w ater flo w , w ater com pressibil-
ity and t her mal expansion which w as t he initial im petus f or t his research.
W it h t he adv ent of lar g e com puting po w er , more recent s tudies of lique-
28
4.1. Back g round
f action phenomena ha v e used discrete element models (DEM) which oper -
ate at t he g r ain scale, and calculate t he equations of mo tion f or t housands
of individuall y tr ack ed cy lindrical or spherical g r ains. In (Goren et al. 2011 )
a 2D DEM model w as coupled t o a continuum model f or fluid flo w , and
t he inter actions of g r ains and fluid w ere calculated f or a sam ple of a f e w
t housand g r ains. Their conditions are typical of 200m t o 2km deep t hin de-
posits sheared at ord( 1) t o ord( 10)m=s
1
which is rele v ant f or f ault g oug e
conditions. Holding boundar y conditions of t heir sam ple at eit her zero fluid
mass flo w , or cons tant fluid pressure, t he y w ere able t o obser v e liquef action
f or bo t h dense and loose sheared assemblies under dr ained and undr ained
conditions. While t heir met hod is in principle applicable t o a wide v ariety of
situations, com puting t he inter actions of lar g e deposits of sand o v er meters
or tens of meters w ould be com putationall y prohibitiv e. Their model does
sho w , ho w e v er , t hat our unders tanding of e v en tablet op sized experiments
ma y be fla w ed, as t he y obser v e liquef action under all conditions in assem-
blies whose bulk density w as eit her relativ el y loose or relativ el y dense. In
an ear lier paper (Goren et al. 2010 ), a nondimensional continuum equation
f or dynamic fluid flo w w as deriv ed which is v alid f or mesoscopic scales and
based on conser v ation of mass t og et her wit h Darcy’ s la w . Their equation
pro vides much of what is required t o anal yze a realis tic soil deposit, t hough
t he y explicitl y neglect certain aspects such as t her mal heating, and t he y do
no t extensiv el y anal yze t he equation in t he context of typical near -surf ace
liquef action conditions.
In t his paper , w e anal yze t he liquef action of satur ated sandy and silty
soils in t he firs t ord( 10) meters belo w t he g round where initial w ater pres-
sures are wit hin a f e w atmospheres of t he t o tal v ertical s tress. Our met hod is
t o deriv e a one dimensional equation f or fluid flo w in t he v ertical direction
based on t he assum p tions of mass conser v ation and Darcy’ s la w f or fluid
flo w t hrough porous media in a manner v er y similar t o Goren et al. ( 2010 ).
Alt hough t he same deriv ation can be triviall y extended t o 3D mat hemati-
call y , it pro vides no additional qualitativ e unders tanding f or t he main point
of t his s tudy , which is t hat fluid flo w and soil inhomog eneity are of critical
im portance in t he liquef action phenomenon. In de v elopment of full pre-
dictiv e models f or realis tic heterog eneous 3D soil deposits ho w e v er , later al
1
N o te on our use of asym p t o tic no tation: x = O( y) means t hatj x= yj< C f or some pos-
itiv e cons tant C which is expected t o be no t extremel y lar g e, x = ord( y) meansj x= yj and
j y= xj should bo t h be treated as close t o 1 and x =o( y) meansj x= yj is negligibl y small.
29
4.1. Back g round
w ater flo w mus t be accounted f or as it will influence t he w ater pressure as
w ell.
By treating a com prehensiv e set of dependencies including tem per ature
effects and t he com pressibility of w ater po tentiall y containing some small
quantity of g as such as arising from or g anic processes or dissol v ed g asses,
w e arriv e at a v er y g ener al result which lea v es no furt her deg rees of freedom
t hat could be expected t o contribute significantl y t o w ater pressurization in
t hese soils. By nondimensionalizing t he model using a dis tinguished limit
(F o w ler 1998 ) w e f ocus attention on t he relativ e size of t he v arious effects
t hat can cause fluid pressurization and tr ansf er of s tress from soil particles
t o t he w ater under t he rele v ant circums tances. By sho wing t hat heating and
com pressibility of w ater are bo t h perturbation effects, our results clarify t hat
g r ains collapsing int o denser arr ang ements pressurizes t he fluid. Our results
also sho w t hat pressure diffusion can no t be neglected in t he liquef action
process. Indeed, w e sho w t hat f or loose sands, diffusion is equall y im portant
t o densification on timescales significantl y shorter t han t he dur ation of an
eart hquak e and e v en shorter t han a single loading cy cle.
F or la y ered deposits or sam ples wit h spatiall y v ar ying per meability , g r a-
dients in t he per meability field are critical in deter mining where liquef action
onset might occur . Though t his has been sho wn in a v ariety of experiments,
it has ne v er been explicitl y anal yzed mat hematicall y t o deter mine ho w im-
portant t he effect is in g ener al. Our model explains t he dominant effects in a
r ang e of conditions from v er y lo w per meability soils typical of silts, t hrough
sands, and e v en in t he case of high per meability g r a v els where diffusion is
more dominant. One v alue of t he model is t hat it predicts certain ne w f ea-
tures, such as lar g e settlements t og et her wit h bulk fluid flo w t o t he surf ace
in loose high per meability deposits, and also sugg es ts a ne w mechanism of
liquef action po tentiall y mediated b y heating wit hin silts e v en if t he silt’ s ini-
tial density is such t hat sus tained contr activ e loading does no t occur . W e
com pute solutions f or se v er al exam ple scenarios t hat mimic t he obser v ed
dynamics in tablet op (K okusho 1999 ; K okusho and K ojima 2002 ) and cen-
trifug e experiments (F ieg el and K utter 1994 )
Based on t he results of t hese anal yses, w e belie v e t hat liquef action e v al-
uation approaches should explicitl y ackno w ledg e t he im portance of fluid
flo w . F uture efforts should be placed on e v aluating per meability , porosity ,
and t he ener g etic po tential f or eart hquak e-induced porosity chang es o v er
entire deposits of soil, no t mainl y t he exis tence of la y ers of loose sands and a
no tional “resis tance t o liquef action” based on shear s tresses, as is common in
30
4.2. Model Cons tr uction
current Geo technical Engineering pr actice (Idriss and Boulang er 2006 ; Seed
et al. 2003 ).
In addition t o a descrip tion of t he fluid flo w , t he prediction of t he g r ain
mo tions is of g reat interes t, as it is a driving f orce f or t he fluid pressuriza-
tion and mig r ation. A t present t he aut hors are una w are of an y continuum
model which could reliabl y subs titute f or t he results of full y coupled DEM
and y et be com putationall y tr actable f or lar g e scales of spatiall y v ariable soil
deposits. Because of t his limitation, w e appl y an assumed g r ain def or ma-
tion, usuall y a cons tant in space and time, and in v es tig ate t he predicted re-
sults in ter ms of t he induced fluid flo w .
4.2 Model Constr uction
W e assume a horizontall y la y ered deposit of soil, and a v ertical coor dinate
axis pointing upw ar d wit h zero ref erence point at a deep im per meable la y er
belo w t he region of interes t. W e choose t his coor dinate sys tem r at her t han
one oriented do wn w ar d from t he g round surf ace since t he g round surf ace
ma y settle do wn w ar d. W e are interes ted in t he firs t f e w tens of meters belo w
t he g round surf ace, as t his is where small pore v olume s tr ains can produce
chang es in w ater pressure t hat can easil y produce a zero effectiv e s tress con-
dition. W e consider a t hin horizontal slice of t his material wit h cross sec-
tional area A which is v er y lar g e com pared wit h t he typical size of a g r ain,
and t hickness dz which is on t he or der of se v er al g r ain diameters. F or t he sat-
ur ated sands w e consider , t he fluid in our differential v olume at pressure P
and density ρ fills t he inter -g r ain spaces whose v olume is V
f
= φ A dz , where
φ is t he bulk porosity (t he fr action of t o tal v olume occupied b y pore space).
Since t he bulk modulus of t he g r ains is typicall y an or der of magnitude
higher t han t hat of t he fluid, which is already quite high, w e assume t hat
t he individual g r ains do no t com press b y an y appreciable amount. Ho w -
e v er , g r ain rearr ang ement can cause chang es in porosity .
While t he usual assum p tion in liquef action research is t hat fluid flo w
does no t occur on t he time scale of t he eart hquak e, w e ins tead model t he
r ate of chang e of fluid mass wit hin t he v olume of a t hin slice b y examining
t he difference in t he fluid flux out of t he t op surf ace and int o t he bo tt om sur -
f ace. W e can es timate t he fluid flo w r ate v
f
relativ e t o t he g r ains mo ving at
v elocity v
g
using Darcy’ s la w v
f
v
g
=
k
d
μ
@ P
@ z
, where k
d
is per meability and
μ is t he fluid viscosity . Darcy’ s la w has been es tablished t o giv e an accur ate
31
4.2. Model Cons tr uction
effectiv e fluid flo w v elocity v
f
which giv es t he proper “apparent” v olumetric
flo w r ate when multiplied b y t he t o t al cross sectional area of a representativ e
v olume. This equation represents an appro ximation of t he full fluid flo w
equation when onl y viscous effects and pressure g r adients dominate.
Assuming t hat t he horizontal slice can be made t hin enough relativ e t o
t he lengt h scale o v er which fluid v elocity v aries, t he resulting differential
equation f or chang e in fluid mass wit hin t he differential slice is:
d
dt
( ρ φ A dz) dt =
@
@ z
(
A ρ
(
k
d
μ
@ P
@ z
+ v
g
)
dt
)
dz Δ M
static
(4.1)
Here t he lef t hand side expresses t he chang e in a short time dt of t he t o tal
fluid mass in t he slice. The right hand side expresses t he difference in fluid
mass flux which is mo ving relativ e t o t he g r ains at v elocity
k
d
μ
@ P
@ z
. T o g et t he
absolute v elocity in a N e wt onian ref erence fr ame w e add v
g
. F or t he actual
v elocity of t he w ater s tream t hrough t he res tricted cross section of v oids w e
w ould multipl y A b y φ and v
w
b y 1= φ on t he right hand side which cancels
and leads t o t he same equation here, but it should be mentioned t hat t he
actual v elocity of t he w ater t hrough t he res tricted cross section is t he one
of interes t f or an y one seeking t o model t he fluid dr ag f or g r ain mo v ement
equations.
Because liquef action of soils tak es place o v er a v er y small r ang e of fluid
density and tem per ature conditions, w e can eliminate t he density ρ from
t he equation b y sol ving f or it in ter ms of pressure and tem per ature using t he
T a y lor expanded s tate equation:
P( ρ; T) = P
0
+
K
B
ρ
0
( ρ ρ
0
)+ K
B
α( T T
0
) (4.2)
Where P
0
, ρ
0
and T
0
are ref erence pressure, density , and tem per ature le v -
els f or t he w ater , K
B
is t he bulk modulus of w ater( V
@ P
@ V T
) , α is t he coefficient
of v olumetric t her mal expansion(
1
V
@ V
@ T P
) , and ρ and T are t he actual density
and tem per ature. The tem per ature dependence can be deriv ed from allo w -
ing a v olume of w ater t o be heated and expand under cons tant pressure, and
t hen be squeezed back t o its original v olume at cons tant tem per ature. F or
our pur poses w e neglect t he v er y small difference betw een t he adiabatic and
iso t her mal bulk modulus of w ater .
The ter m Δ M
static
is a correction necessar y t o account f or t he f act t hat
t here is no chang e in t he slice’ s w ater mass when s tatic g r a vitational condi-
32
4.2. Model Cons tr uction
tions appl y . F or exam ple, due t o g r a vity t here is a s tatic pressure g r adient.
Also, t he per meability ma y chang e wit h dep t h, and w ater density chang es
minutel y wit h dep t h. Setting t he lef t side of t he equation equal t o zero when
s tatic g r a vitational conditions appl y and sol ving f or t he correction giv es
Δ M
static
=
ρ
2
g A
μ
(
g k
d
ρ
0
K
B
k
d
ρ
dρ
dz
@ k
d
@ z
)
dt dz (4.3)
In equations ( 4.1 ) and ( 4.3 ), f or a v ertical cartesian axis z , A , dz , and dt
are cons tant, k
d
depends on location z , μ depends on tem per ature T , ρ is
deter mined from P and T using t he s tate equation, and P , T and φ depend
on time t and position z .
N o te t hat w e ma y com pensate f or t he presence of a small fr action of t he
v oids containing g as b y adjus ting K
B
. If a v olume of w ater contains small
bubbles occup ying a small fr action ε of t he v olume, t hen w e can appro ximate
t he com pressibility ( 1= K
B
) using V
tot
= V
W
+ V
g
and 1= K
Bmix
=
1
V
tot
@ V
tot
@ P
t og et her wit h t he s tate equation V
g
= K T= P f or an ideal g as and t he T a y lor
series s tate equation f or t he liquid V
w
= V
0
( 1 ( P P
0
)= K
B
) . Assuming
iso t her mal g a s com pression due t o t he t her mal equilibrium wit h t he w ater ,
V
g
= V
0
P
0
= P . In a mixture wit h ε fr action of g as, t his leads t o t he equation:
V
tot
( P) = ( 1 ε) V
tot 0
( 1
( P P
0
)
K
B
)+ ε V
tot 0
P
0
P
e v aluating
K
Bmix
=
1
V
tot 0
@ V
tot
@ P
at P
0
, w e g et:
P
0
K
B
mix
= ε
(
1
P
0
K
B
)
+
P
0
K
B
This is a s trong function of ε since P
0
= K
B
is small f or atmospheric scale pres-
sures. This ma y be rele v ant in situations in v ol ving g as bubbles f or med b y
biological or g anisms or recent r ains bearing an ex cess of dissol v ed g as, or in
tidal areas where flux es of w ater and g as bubbles are lar g e.
Since w e used t he s tate equation in v ol ving P and T t o eliminate ρ , our
equation no w includes t he tem per ature T and in expanded f or m will include
bo t h
@ T
@ t
and
@ T
@ z
. W e use t he f ollo wing heat equation f or a flo wing fluid t o
relate t hese quantities:
33
4.2. Model Cons tr uction
dT
dt
+
@ T
@ z
v
w
=
1
ρ
avg
c
v
(
k
T
@
2
T
@ z
2
+ T
@ S
@ t
Q
)
(4.4)
In ( 4.4 ) c
v
is a bulk a v er ag ed specific heat capacity , ρ
avg
is t he bulk a v er ag e
density , k
T
is t he effectiv e t her mal conductivity , and
@ S
@ t Q
is t he chang e in
entrop y density due t o heat g ener ation.
In addition, t he conser v ation of g r ain mass w as used t o replace
@ v
g
@ z
in
equation 4.1 wit h a ter m in v ol ving
@ φ
@ t
and a con v ectiv e ter m in v ol ving v
g
,
which will become a perturbation due t o t he small settlement v elocities.
@( 1 φ)
@ t
=
@( 1 φ) v
g
@ z
!
@ v
g
@ z
=
1
( 1 φ)
(
@ φ
@ z
v
g
+
@ φ
@ t
)
(4.5)
In t his f or m, equation 4.1 is no t y et useful f or anal ysis and calculation,
as wit hin each deriv ativ e t here are multiple v ariables whose v alues are de-
pendent on time or position. Our next task will be t o expand t he fluid equa-
tion ( 4.1 ) t hrough t he com prehensiv e application of t he chain r ule, so w e
can arriv e at an equation f or t he r ate of chang e of pressure
@ P
@ t
and t hen mak e
t his equation dimensionless t o help e v aluate t he relativ e im portance of v ar -
ious ter ms. Alt hough it is tem p ting at t his point t o tr y t o sim plify t his equa-
tion b y making assum p tions about t he negligibility o f certain effects such as
com pressibility of w ater , t he heat g ener ation r ate, and t he lik e t o reduce t he
com plexity of t he expanded equation, w e belie v e t hat such ear l y sim plifying
assum p tions ha v e led t he field as tr a y in g ener al, and w e ins tead choose t o
expand t he equation full y and deal wit h t he negligible ter ms onl y af ter a full
dimensionless anal ysis.
W e em plo y ed t he freel y a v ailable Maxima com puter alg ebr a sys tem (The
Maxima Project 2011 ) t o calculate t he expanded v ersion of equation 4.1 . By
identifying and dropping negligible ter ms w e arriv e at a sim plified “full”
equation which cap tures all of t he dynamics of interes t including firs t or -
der perturbations (equation 4.18 is t he full y anal yzed dimensionless v ersion
belo w).
In or der t o arriv e at t he sim plified equation, w e needed t o unders tand t he
dominant balances of effects a nd identify t he negligible ter ms. T o do t his a
proper nondimensionalization is required. The choice of t he char acteris tic
scales f or nondimensionalization is one of t he critical aspects of our anal ysis,
since it f ocuses t he equation on t he specific situations of interes t.
34
4.2. Model Cons tr uction
T o f ocus our unders tanding on g eo technicall y rele v ant scales, w e chose
scales f or P , k
d
and φ b y ref erence t o t he situation in t he firs t se v er al meters
of sandy soil, where P= P
0
= O( 1) t hroughout. W e introduce k
d0
which rep-
resents a typical or der of magnitude f or per meability of moder atel y coarse
loose sands, such t hat t he per meability r ang e f or mos t sands f alls wit hin
about a f act or of 4 of t his v alue. W e also introduce φ
0
which represents a
typical v alue f or porosity in a loose sand. Gener al sands are betw een about
70% and 110% of t his φ v alue. These v alues w ere chosen b y ref erence t o
tables giv en in Bear d and W e y l ( 1973 ).
P
0
= P
atm
= 100kPa
k
d0
= 1: 5 10
10
m
2
φ
0
= 0: 44
(4.6)
W e wro te equations f or t he dimensional v ariables in ter ms of nondimen-
sional “primed” or “hatted” v ariables as f ollo w s:
φ = φ
0
(
1+
Δ φ
φ
0
φ
′
)
z = Δ z z
′
v
g
=
6Δ φ Δ z
Δ t
v
′
g
^ φ = φ= φ
0
k
d
= k
d0
k
′
d
S = Δ S S
′
P = P
0
+ P
0
P
′
t = Δ t t
′
T = T
0
(
1+
Δ S
C
v0
T
′
)
(4.7)
Here C
v0
1594J=kg=K represents a mass w eighted a v er ag e of t he specific
heat capacity of g r ains and w ater . W e t ook numerical v alues f or t he heat ca-
pacities of w ater and g r ains as C
vw
= 4156: 7J=kg=K and C
vgr
= 835J=kg=K .
W e t hen subs tituted t hese expressions int o t he expanded dimensional f or m
of equation 4.1 and sol v ed f or t he dimensionless
@ P
′
@ t
′
.
During liquef action, due t o t he high bulk modulus of w ater , φ does no t
need t o chang e v er y much in or der t o induce a unit of dimensionless pres-
sure. W e ha v e introduced t he v ariable φ
′
which chang esO( 1) when φ chang es
b y t his tin y amount required f or pressurization. On t he o t her hand ^ φ repre-
sents a nor malized v alue of φ as a fr action of its typical size, a v alue which
can be considered v er y near l y cons tant during t he e v ent. The scaling of di-
mensional v elocity v
g
w as chosen so t hat 1 unit of v
′
g
represents a g round
settlement of about 10 cm during a 30 s eart hquak e af ter com puting Δ φ as
described belo w .
The remaining scale cons tants Δ z , Δ t , Δ S and Δ φ can be chosen arbitr aril y ,
ho w e v er t he situation wit h t he maximal inter pla y of different effects is t he
35
4.2. Model Cons tr uction
so called dis tinguished limit (F o w ler 1998 ), where all of t he rele v ant ter ms in
t he equation ha v e coefficient equal t o 1 under typical conditions so t hat one
unit of chang e in an y v ariable or deriv ativ e induces t he same size of effect
com pared t o chang es of an y o t her v ariable or deriv ativ e. By calculating our
scale f act ors t o f orce t he equation int o t his f or m, w e identify t he regime in
which maximal inter pla y occurs.
T o calculate t he scale cons tants, w e e v aluate t he coefficients under t he
conditions k
′
d
= 1 and φ = φ
0
as representativ e of typical conditions, and
t hen require t he im portant ter ms
@
2
P
′
@ z
′ 2
,
@ φ
′
@ t
′
, and
@ P
′
@ z
′
@ k
′
d
@ z
′
which represent t he
effect of pressure diffusion, g r ain sk elet on contr action, and spatial v ariability
in per meability respectiv el y t o ha v e coefficients equal t o 1. In addition w e
require t he initial dimensionless g r a vitational s tatic pressure g r adient t o be
1. These f our equations are sufficient t o cons tr ain t he v alues of Δ z , Δ t , Δ S and
Δ φ and w e sol v e f or t hem as f ollo w s:
Δ z = P
0
= ρ
0
g 10m
Δ t =
φ
0
μ P
2
0
g
2
k
d0
K
B
ρ
2
0
0: 14s
Δ S =
C
v0
P
0
γ
S
α K
B
T
0
6: 0 10
3
J=kg=K( γ
S
= 1)
Δ φ =
φ
0
( 1 φ
0
)
( 2 φ
0
)
P
0
K
B
7: 2 10
6
(4.8)
Clear l y our assum p tion t hat diffusion, contr action, and per meability v ari-
ation are all im portant in realis tic flo w regimes is v alidated b y t he magnitude
of t hese scale cons tants. The v alue of Δ t , which is a diffusion timescale re-
lated t o( Δ z)
2
= D where D is t he pressure diffusivity , is small com pared t o an
eart hquak e dur ation of around 30 s, and e v en t o a single cy cle of loading of
appro ximatel y 1 s. The lengt h scale Δ z is O( 1) relativ e t o t he size of deposits
of interes t, and Δ φ is a reasonabl y achie v able fluctuation in φ . The f act t hat
Δ t is small com pared t o t he eart hquak e dur ation or t he single loading cy cle
means t hat a near -s teady -s tate appro ximation of t he flo w should be v alid f or
predicting t he pressure field at t he end of eart hquak e shaking. This is tr ue
e v en t hough t he f orcing ma y chang e significantl y o v er a loading cy cle as t he
diffusion ter m
@
2
P
′
@ z
′ 2
im plies a g aussian w eighted a v er aging o v er one dimen-
sionless unit of lengt h during one dimensionless unit of time. Theref ore all
points in t he deposit are influencing all o t her points and no localized pres-
36
4.2. Model Cons tr uction
sure chang es can las t more t han O( 0: 14)s wit hout bo t h r apid and sus tained
g r ain sk elet on contr action in a single direction t o balance t he diffusion. W e
expect t hat dr ag f orces on particles pre v ent t he sus tained contr action in a sin-
gle direction from becoming t oo r apid (as will be sho wn in later exam ples),
and since t he typical eart hquak e loading fluctuates in bo t h directions, t he
f as t fluctuations in loading are a v er ag ed out and t he pressure field remains
near s teady s tate when time a v er ag ed o v er at leas t a cy cle.
N o te in t he abo v e scales w e ha v e lef t γ
S
as a nondimensional par ame-
ter which measures t he fr action of t he “dis tinguished limit” le v el of heating
which actuall y occurs in pr actice. If γ
S
= 1 , t hen t his im plies a char acteris tic
tem per ature chang e r ate of 1: 05K=s which f or a 30s eart hquak e im plies 32K
tem per ature chang e, a clear l y untenable quantity : t he sand w ould be scald-
ing ho t, a phenomenon no t reported pre viousl y . The dis tinguished limit
wit h heating is t heref ore more rele v ant f or some mechanical appar atus in
which exter nal heating is applied, or per haps a g eo t her mal or f ault goug e
condition or wit hin r apidl y shearing bands where lar g e quantities of ener gy
are being dissipated in small spatial regions, such as in landslides or t hin
partiall y liquefied la y ers. Ho w e v er , as pointed o ut b y N emat-N asser and
Shok ooh ( 1979 ) and v erified in field obser v ations b y Da vis and Berrill ( 2001 ),
dissipated mechanical ener gy is highl y correlated wit h w ater pressure de v el-
opment. Clear l y , t he w a v e ener gy is used b y t he g r ains as activ ation ener gy
f or rearr ang ement, and t his process is irre v ersible and mus t produce some
heat. T o es timate t he actual scale of heating wit hin t he soil during t he ear l y
s tag es of typical liquef action e v ents, and t hereb y cons tr ain γ
S
, w e made ref-
erence t o t he entrop y g ener ated b y t he flo w of w ater t hrough a char acteris-
tic soil under a char acteris tic pressure g r adient. W e calculate t his entrop y
scale b y assuming t hat t he w ater flo w s horizontall y f or Δ t seconds t hrough
a cons tant pressure g r adient P
0
= Δ z at a v elocity giv en b y Darcy’ s la w , and
t he w or k done b y t he pressure field (t h e chang e in ent halp y) is com pletel y
con v erted t o entrop y at tem per ature T
0
, via dr ag from t he particles:
γ
S
min
=
(
P
0
Δ z
)(
k
d0
μ
P
0
Δ z
)
Δ t
T
0
ρ
w
(
1
Δ S
γ
S
= 1
)
=
α φ
0
P
atm
C
v0
ρ
w
6 10
6
(4.9)
This is a lo w er bound f or t he relativ e entrop y g ener ation in t he soil as it
37
4.2. Model Cons tr uction
does no t include an y effect of t he g r ain mo tion. T o bound t he r ang e of rea-
sonable v alues f or γ
S
w e assumed t hat during a 30 second eart hquak e soil
tem per ature should no t chang e more t han ord( 1)K since an or der of mag-
nitude more w ould no t be ignor able, and tem per ature mus t chang e at leas t
ord( 10
4
)K due t o t he flo w entrop y calculated abo v e. As a representativ e
v alue, w e t heref ore chose γ
S
5 10
3
under t he assum p tion t hat t he g r ain
mo tion contributes significantl y t o entrop y production but t hat t he tem per a-
ture chang e is s till a fr action of a K el vin. The im plied char acteris tic chang e in
tem per ature f or t he soil o v er all is 0: 24K in 30s . This w ould seem t o be a lar g e
fr action of t he amount t hat could occur while s till being neglected in mos t
field or labor at or y experiments. Ev en if t he tr ue v alue is an or der of magni-
tude lar g er t han t his, our treatment of t he effect as a small perturbation will
remain v alid.
Af ter nondimensionalizing, w e replaced each of f our im portant nondi-
mensional g roups wit h a single nondimensional symbol f or each g roup as
f ollo w s:
γ
vol
=
P
atm
K
B
(4.10)
γ
μ
T
=
1
α μ
d μ
dT
(4.11)
G
gr
= ρ
grain
= ρ
0
(4.12)
β(^ φ) =
1+
1
( 1 φ
0
^ φ)
^ φ
(
1+
1
( 1 φ
0
)
)=
( φ
0
1)( φ 2)
( φ
0
2)^ φ( φ 1)
(4.13)
A bo v e, γ
vol
expresses t he com pressibility of w ater in ter ms of t he v olumet-
ric s tr ain required t o g ener ate a char acteris tic unit of pressure, γ
μT
repre-
sents fr actional chang e in viscosity f or a nondimensional unit of tem per ature
chang e, and G
gr
is t he relativ e density of g r ains t o w ater , typicall y about 2: 6
t o 2: 7 ; here, w e ha v e used 2: 655 . β is t he local nondimensional bulk modulus
of t he mixture. It giv es t he dimensionless pressure induced b y a unit chang e
of φ
′
in t he absence of flo w , and is a function of t he local ^ φ v alue r anging in
pr actice betw een about 1 and 2 f or realis tic v alues of ^ φ 2 [ 0: 5; 1: 1] . F or im-
portant cons tants related t o w ater , w e used t he v alues of density , heat capac-
ity , viscosity , and o t her im portant properties a v ailable online from t he NIS T
database of t he properties of w ater f or a ref erence pressure of P=100kP a and
tem per ature T=293.15K ( 20
◦
C )(NIS T 2012 ).
38
4.2. Model Cons tr uction
By cons tr uction, f or typical conditions of soil liquef action all t he v ari-
ables and t heir deriv ativ es are O( 1) so t hat t he relativ e size of each ter m is
deter mined from its coefficient. In certain conditions t hese assum p tions are
violated, such as near sudden tr ansitions in per meability . T o eliminate t he
negligible ter ms and identify t he leading or der ter ms of t he perturbation w e
perf or med an asym p t o tic anal ysis of t he equation. R etained perturbation
ter ms are rele v ant onl y near lar g e g r adients in per meability , lar g e g r adients
in pressure, or extremel y lo w per meability regions.
T o see t he de v elopment of t he full y anal yzed equation w e firs t sho w t he
full expanded dimensional equation:
39
4.2. Model Cons tr uction
dP
dt
=
@ k
d
@ z
(
K
B
(
α
2
g ρ
0
( T
0
T)
2
μ φ
+
2 α g ρ
0
( T
0
T)
μ φ
+
g ρ
0
μ φ
)
+ α
(
2 g ρ
0
P( T
0
T)
μ φ
2 g P
atm
ρ
0
( T
0
T)
μ φ
)
+
g ρ
0
P
2
μ φ
2 g P
atm
ρ
0
P
μ φ
+
g P
2
atm
ρ
0
μ φ
μ φ
+
2 g ρ
0
P
μ φ
2 g P
atm
ρ
0
μ φ
)
α
2
g
2
k
d
ρ
2
0
( T
0
T)
2
μ φ
+
@ P
@ z
@ T
@ z
(
K
B
(
α
(
k
d
( φ 1)
μ φ
k
d
d μ
dT
( T
0
T)
μ
2
φ
)
k
d
d μ
dT
μ
2
φ
)
+
d μ
dT
(
k
d
P
atm
μ
2
φ
k
d
P
μ
2
φ
)
)
+
α
(
2 g
2
k
d
P
atm
ρ
2
0
( T
0
T)
μ φ
2 g
2
k
d
ρ
2
0
P( T
0
T)
μ φ
)
2 g
2
k
d
ρ
2
0
P
μ φ
+
2 g
2
k
d
P
atm
ρ
2
0
μ φ
μ φ
+
@ T
@ z
(
K
B
(
α
(
μ v
g
g k
d
ρ
0
)
μ φ
α
2
g k
d
ρ
0
( T
0
T)
μ φ
)
+ α
(
g k
d
P
atm
ρ
0
μ φ
g k
d
ρ
0
P
μ φ
)
)
+
@ φ
@ t
(
K
B
(
α ( φ 2) ( T
0
T)
( φ 1) φ
φ 2
( φ 1) φ
)
( φ 2) P
( φ 1) φ
+
( φ 2) P
atm
( φ 1) φ
)
+
@
2
P
@ z
2
(
K
B
(
α k
d
( T
0
T)
μ φ
+
k
d
μ φ
)
+
k
d
P
μ φ
k
d
P
atm
μ φ
)
+
@ k
d
@ z
@ P
@ z
(
K
B
(
α ( T
0
T)
μ φ
+
1
μ φ
)
+
P
μ φ
P
atm
μ φ
)
+ K
B
0
@
α
(
C
v 0
@ φ
@ z
C
′
v
(^ φ) v
g
T
0
+ φ
2 @ S
@ t
T φ
@ S
@ t
T C
v 0
@ φ
@ z
C
′
v
(^ φ) v
g
T
)
C
v 0
( φ 1) φ C
′
v
(^ φ)
+
@ φ
@ z
v
g
( φ 1) φ
1
A
+
@ P
@ z
(
α g k
d
ρ
0
( T
0
T)
μ φ
+
g k
d
ρ
0
P
μ φ
g k
d
P
atm
ρ
0
μ φ
μ φ
μ v
g
g k
d
ρ
0
μ φ
)
2 α g
2
k
d
ρ
2
0
( T
0
T)
μ φ
+
α K
T
K
B
@
2
T
@ z
2
C
v 0
C
′
v
(^ φ)
(
φ ρ
0
φ ρ
gr
+ ρ
gr
)+
k
d
(
@ P
@ z
) 2
μ φ
+
g
2
k
d
ρ
2
0
P
2
μ φ
+
2 g
2
k
d
P
atm
ρ
2
0
P
μ φ
g
2
k
d
P
2
atm
ρ
2
0
μ φ
μ φ
+
μ
@ φ
@ z
v
g
P g
2
k
d
φ ρ
2
0
+ g
2
k
d
ρ
2
0
μ ( φ 1) φ
@ φ
@ z
P
atm
v
g
( φ 1) φ
(4.14)
Af ter replacing t he dimensionless g roups wit h t heir dimensionless pa-
r ameters, and reor g anizing, t he full dimensionless equation is:
40
4.2. Model Cons tr uction
dP
′
dt
′
=
@ S
′
@ t
′
T
′
γ
2
S
γ
vol
α C
′
v
(^ φ) T
0
+ γ
S
γ
2
vol
0
@
k
′
d
φ
0
@ P
′
@ z
′
(
P
′ @ T
′
@ z
′
γ
μ
T
+ T
′
)
φ
6
(
φ
0
1
)
2
φ
2
0
@ φ
′
@ z
′
T
′
v
′
g
( φ 1) φ
(
φ
0
2
)
2
k
′
d
φ
0
(
P
′ @ T
′
@ z
′
2 T
′
)
φ
2
@ k
′
d
@ z
′
φ
0
P
′
T
′
φ
1
A
+ γ
S
γ
vol
0
@
k
′
d
φ
0
@ P
′
@ z
′
@ T
′
@ z
′
(
γ
μ
T
φ+ 1
)
φ
+
6
(
φ
0
1
)
φ
0
@ T
′
@ z
′
v
′
g
φ
(
φ
0
2
)
k
′
d
φ
0
@ T
′
@ z
′
φ
k
′
d
φ
0
@
2
P
′
@ z
′
2
T
′
φ
@ k
′
d
@ z
′
φ
0
@ P
′
@ z
′
T
′
φ
+
( φ 2)
(
φ
0
1
)
φ
0
@ φ
′
@ t
′
T
′
( φ 1) φ
(
φ
0
2
)
2
@ k
′
d
@ z
′
φ
0
T
′
φ
1
A
+ γ
2
S
γ
2
vol
(
k
′
d
φ
0
@ P
′
@ z
′
T
′ @ T
′
@ z
′
γ
μ
T
φ
+
k
′
d
φ
0
T
′ @ T
′
@ z
′
φ
+
@ k
′
d
@ z
′
φ
0
T
′
2
φ
)
k
′
d
φ
0
T
′
2
γ
2
S
γ
3
vol
φ
+
2 k
′
d
φ
0
P
′
T
′
γ
S
γ
3
vol
φ
k
′
d
φ
0
P
′
2
γ
3
vol
φ
+
(
6
(
φ
0
1
)
2
φ
2
0
@ φ
′
@ z
′
P
′
v
′
g
( φ 1) φ
(
φ
0
2
)
2
+
k
′
d
φ
0
P
′ @ P
′
@ z
′
φ
+
@ k
′
d
@ z
′
φ
0
P
′
2
φ
2 k
′
d
φ
0
P
′
φ
)
γ
2
vol
+
0
@
6
(
φ
0
1
)
φ
0
@ P
′
@ z
′
v
′
g
φ
(
φ
0
2
) +
6
(
φ
0
1
)
2
φ
2
0
@ φ
′
@ z
′
v
′
g
( φ 1) φ
(
φ
0
2
)
2
+
k
′
d
φ
0
P
′ @
2
P
′
@ z
′
2
φ
+
k
′
d
φ
0
(
@ P
′
@ z
′
) 2
φ
+
@ k
′
d
@ z
′
φ
0
P
′ @ P
′
@ z
′
φ
+
k
′
d
φ
0
@ P
′
@ z
′
φ
( φ 2)
(
φ
0
1
)
φ
0
@ φ
′
@ t
′
P
′
( φ 1) φ
(
φ
0
2
) +
2
@ k
′
d
@ z
′
φ
0
P
′
φ
k
′
d
φ
0
φ
1
A
γ
vol
γ
thermC
@
2
T
′
@ z
′
2
γ
S
(
G
gr
φ φ G
gr
)
C
′
v
(^ φ)
+
@ S
′
@ t
′
γ
S
C
′
v
(^ φ)
+
k
′
d
φ
0
@
2
P
′
@ z
′
2
φ
+
@ k
′
d
@ z
′
φ
0
@ P
′
@ z
′
φ
( φ 2)
(
φ
0
1
)
φ
0
@ φ
′
@ t
′
( φ 1) φ
(
φ
0
2
) +
@ k
′
d
@ z
′
φ
0
φ
(4.15)
Sim plifying t his equation b y taking t he firs t or der T a y lor series in t he
small par ameters produces:
41
4.2. Model Cons tr uction
dP
′
dt
′
=
0
@
6
(
φ
0
1
)
φ
0
@ P
′
@ z
′
v
′
g
φ
(
φ
0
2
) +
6
(
φ
0
1
)
2
φ
2
0
@ φ
′
@ z
′
v
′
g
( φ 1) φ
(
φ
0
2
)
2
+
k
′
d
φ
0
P
′ @
2
P
′
@ z
′
2
φ
+
k
′
d
φ
0
(
@ P
′
@ z
′
) 2
φ
+
@ k
′
d
@ z
′
φ
0
P
′ @ P
′
@ z
′
φ
+
k
′
d
φ
0
@ P
′
@ z
′
φ
( φ 2)
(
φ
0
1
)
φ
0
@ φ
′
@ t
′
P
′
( φ 1) φ
(
φ
0
2
)
+
2
@ k
′
d
@ z
′
φ
0
P
′
φ
k
′
d
φ
0
φ
1
A
γ
vol
+
@ S
′
@ t
′
γ
S
C
′
v
(^ φ)
+
k
′
d
φ
0
@
2
P
′
@ z
′
2
φ
+
@ k
′
d
@ z
′
φ
0
@ P
′
@ z
′
φ
( φ 2)
(
φ
0
1
)
φ
0
@ φ
′
@ t
′
( φ 1) φ
(
φ
0
2
) +
@ k
′
d
@ z
′
φ
0
φ
(4.16)
Due t o t he lar g e size of γ
μT
t he combination γ
S
γ
μT
can no t be consid-
ered second or der in size, and w e res t ore ter ms related t o t his combination
yielding:
dP
′
dt
′
= γ
vol
0
@
k
′
d
φ
0
@ P
′
@ z
′
@ T
′
@ z
′
γ
S
γ
μ
T
φ
6
(
φ
0
1
)
φ
0
@ P
′
@ z
′
v
′
g
φ
(
φ
0
2
)
+
6
(
φ
0
1
)
2
φ
2
0
@ φ
′
@ z
′
v
′
g
( φ 1) φ
(
φ
0
2
)
2
+
k
′
d
φ
0
P
′ @
2
P
′
@ z
′
2
φ
+
k
′
d
φ
0
(
@ P
′
@ z
′
) 2
φ
+
@ k
′
d
@ z
′
φ
0
P
′ @ P
′
@ z
′
φ
+
k
′
d
φ
0
@ P
′
@ z
′
φ
( φ 2)
(
φ
0
1
)
φ
0
@ φ
′
@ t
′
P
′
( φ 1) φ
(
φ
0
2
) +
2
@ k
′
d
@ z
′
φ
0
P
′
φ
k
′
d
φ
0
φ
1
A
+
@ S
′
@ t
′
γ
S
C
′
v
(^ φ)
+
k
′
d
φ
0
@
2
P
′
@ z
′
2
φ
+
@ k
′
d
@ z
′
φ
0
@ P
′
@ z
′
φ
( φ 2)
(
φ
0
1
)
φ
0
@ φ
′
@ t
′
( φ 1) φ
(
φ
0
2
) +
@ k
′
d
@ z
′
φ
0
φ
(4.17)
The final s tag es of anal ysis are an asym p t o tic anal ysis of t he perturbation
ter ms (t hose multiplied b y γ
vol
) t o retain onl y t hose ter ms which are leading
or der . In or der f or t he perturbations t o be im portant, w e mus t ha v e t he un-
perturbed equation ( 4.19 ) in near -equilibrium (o t her wise it will dominate).
Since all t he v ariables t hemsel v es are O( 1) e v er ywhere, in or der f or t he per -
turbation ter ms t o be appreciable, w e need one of t he deriv ativ es t o be lar g e.
F or exam ple, if
@ k d
′
@ z
′
=
1
ε
is lar g e w e can zoom int o t his small spatial region b y
42
4.2. Model Cons tr uction
rescaling z
′
! ε^ z . The dominant ter ms wit h t he right or der of magnitude t o
balance are t hose wit h a f act or of 1= ε
2
under t his rescaling which includes
t he
@
2
P
′
@ z
′ 2
, (
@ P
′
@ z
′
)
2
, and
@ k
′
d
@ z
′
@ P
′
@ z
′
ter ms, as w ell as t he pre viousl y mentioned ter m
in v ol ving γ
μT
, hence t he y are retained. All remaining perturbation ter ms are
dropped, and t his, t og et her wit h some reor g anization of ter ms and sim pli-
fying subs titutions leads t o our final equation:
dP
′
dt
′
= β(^ φ)
@ φ
′
@ t
′
+
1
^ φ
(
k
′
d
@
2
P
′
@ z
′
2
+
@ k
′
d
@ z
′
(
@ P
′
@ z
′
+ 1
))
+
@ S
′
@ t
′
γ
S
C
′
v
(^ φ)
+ γ
vol
(
β(^ φ)
@ φ
′
@ t
′
P
′
+
1
^ φ
(
k
′
d
@ P
′
@ z
′
@ T
′
@ z
′
γ
S
γ
μ
T
+ k
′
d
P
′
@
2
P
′
@ z
′
2
+ k
′
d
(
@ P
′
@ z
′
)
2
+
@ k
′
d
@ z
′
P
′
@ P
′
@ z
′
))
(4.18)
Equation ( 4.18 ) describes all rele v ant f eatures of our model including
small firs t or der perturbations due t o w ater com pressibility and t her mal ef-
f ects. In t his equation C
′
v
(^ φ) is t he nondimensional relativ e heat capacity as
a fr action of C
v0
which chang es due t o spatial chang es in porosity , and β(^ φ)
is t he nondimensional bulk modulus of t he mixture defined abo v e (equation
4.13 ). In equation ( 4.18 ) t he eliminated ter ms all had coefficients smaller t han
about 10
9
.
It is at t his point t hat t he assum p tion of “fluid incom pressibility , ” which
is of ten in v ok ed from t he beginning in soil mechanics, can be applied in a
r ational manner b y taking a limit as γ
vol
! 0 . This assum p tion is a helpful
one and jus tified f or mos t pur poses in v ol ving soil liquef action. The primar y
conditions in which com pressibility can become significant are high pressure
g r adients in lo w per meability soils, such as at t he center of t hin la y ers of
silt, t he spatial regions where t hese conditions occur are necessaril y small
as discussed in t he oscillating per meability exam ple belo w . Thus t he ter m
multiplied b y γ
vol
is dropped, and t he resulting sim plified equation is:
(4.19)
dP
′
dt
′
= β(^ φ)
@ φ
′
@ t
′
+
1
^ φ
(
k
′
d
@
2
P
′
@ z
′
2
+
@ k
′
d
@ z
′
(
@ P
′
@ z
′
+ 1
))
+
@ S
′
@ t
′
γ
S
C
′
v
(^ φ)
43
4.3. Qualitativ e R esults
4.3 Qualit ativ e R esults
Some basic sim plified models f or t he qualitativ e beha vior expected from
t hese equations are no w considered. Alt hough researchers such as Zienkie wicz,
Chang, and Bettess ( 1980 ) ha v e pre viousl y in v es tig ated t he role of dr ainag e
and sho wn t hat dr ainag e is usuall y im portant, t hese results ha v e no t g ener -
all y been appreciated. Here w e examine se v er al k e y experiments perf or med
b y o t hers in t he context of t he predictions from our equations. F irs t, recall
t hat t he char acteris tic time scale Δ t 0: 14s f or t his sys tem is smaller t han
t he time scale f or a typical eart hquak e loading period of about 1s and much
smaller t han a typical eart hquak e dur ation of around 10 60s . This means
t hat t he beha vior of t he w ater pressure under realis tic loading is alw a ys
near l y in s teady s tate. The diffusion ter m in t he equation im plies g aussian
a v er aging of t he pressure field o v er t he entire deposit in time Δ t = 0: 14s so
t hat an y localized pressure f eature will spread o v er t he whole deposit in t his
short time via pressure diffusion induced b y imbalances in flo w described b y
t he right hand side of equation 4.1 . An yord( 1) imbalance in an y of t he ter ms
of equation 4.19 w ould induceord( 1) bars of pressure chang e o v er 0: 14s , re-
tur ning t he sys tem r apidl y t o near s teady s tate. Ev en f or fine sands where
k
d
k
d0
= 20 t he time scale w ould be appro ximatel y t he same or der as a sin-
gle loading cy cle and significantl y shorter t han typical eart hquak e dur ation.
Onl y when per meability is similar t o t he per meability of silts k
d0
= 1000 can
w e assume t hat t he pressure diffusion time is long com pared t o eart hquak e
dur ation. U nder t hose conditions t he pressure chang es are deter mined b y
t he tin y v olumetric chang es caused b y squeezing t he near l y incom pressible
w ater , and t he “undr ained” appro ximation
@ P
′
@ t
′
β(^ φ)
@ φ
′
@ t
′
might be reason-
able. Ev en in regions of silty soils, as w e shall see in exam ples, when a silt
la y er is small and in contact wit h a sandy la y er , v er y high cur v ature of t he
pressure field and lar g e g r adients g ener ated near t hese interf aces can cause
t he silt t o reach s teady s tate pressures on t he same time scale as t he sand.
Clear l y flo w can no t be neglected in eart hquak e liquef action of sand. In f act
f or sands quite t he opposite case occurs, t he flo w required t o induce s teady
s tate pressures happens extremel y quickl y .
T o support t his anal ysis wit h ph ysical data, w e can ref erence f or exam-
ple t he centrifug e model B of F ieg el and K utter ( 1994 ) in v ol ving a fine sand
o v er lain b y silt. Shaking s tarts at appro ximatel y 2 seconds, and las ts until
appro ximatel y 12 seconds. By 4 seconds t he pore pressure in t he sand has
come t o a locall y time-a v er ag ed s teady s tate wit h a f as t oscillation superim-
44
4.3. Qualitativ e R esults
posed due t o continued shaking, and t he silt has a similar if slightl y slo w er
response, achie ving a locall y time-a v er ag ed s teady s tate at appro ximatel y
10 seconds which is bef ore t he end of shaking at 12 seconds (times giv en
in pro t o type dimension). A s teady s tate solution of equation 4.19 based on
time-a v er ag ed r ate of chang e of porosity could be expected t o be quite accu-
r ate at 10s, e v en bef ore t he shaking s t opped, if t he loading w ere tak en t o be
t he time a v er ag e of t he actual loading o v er a f e w seconds.
Because of t his near -s teady s tate, w e can appro ximate t he pressure field
at liquef action b y sol ving f or t he s teady s tate of equation ( 4.19 ) numericall y
t o find t he a v er ag e f orcing
@ φ
′
@ t
′
required t o mak e t he w ater pressure field tan-
g ent t o t he t o tal v ertical s tress. Once w e ha v e t his solution w e can find t he
location of tang ency which is where t he onset of liquef action occurs. W e use
@ φ
′
@ t
′
cons tant in space ex cep t where no ted belo w . Since heating w as ear lier
deter mined t o be a small perturbation f or sands w e ignore it here f or sim-
plicity ex cep t in t he exam ple in v ol ving t her mal liquef action in silts. Also,
w e set β(^ φ) = ^ φ = 1 and ignore t he small effect of spatial v ariability in φ ex-
cep t as it influences k
d
which is spatiall y v ar ying as sho wn in each figure. By
adjus ting t he char acteris tics of t he soil, especiall y t he per meability field and
t he spatial v ariation in t he f orcing, w e can deter mine ho w t he per meability
and per meability g r adient affect t he onset of liquef action.
Exam ple 1: U nif or m Loose Hydraulic Fill A model of a unif or m loose h y -
dr aulic fill such as t hose f ound at reclaimed land in ba ys and port f acilities
w or ldwide w ould be a unif or m deposit wit h k
′
d
= 1 of dep t h z
′
1 . W e
assume f or sim plicity t hat t he w ater table is at t he surf ace of t he g round,
and w e appl y a 0.12 nondimensional pressure unit load on t he g round de-
signed t o model 1/2 meter of concrete pad t o ensure t hat t he effectiv e s tress
is no t zero at t he surf ace (figure 4.1 ). In t his situation, due t o t he unif or m
per meability , t he ter m in equation 4.19 in v ol ving t he s patial deriv ativ e of k
d
disappears. If w e assume t he contr activ e f orcing
@ φ
′
@ t
′
is cons tant in space and
w e are sol ving f or s teady s tate conditions, t hen t he cur v ature of t he pressure
field is cons tant and w e ha v e a par abolic pressure cur v e.
In figure 4.1 w e sho w t he results of t his simulation and w e find t hat t he
soil jus t initiates liquef action via loss of g r ain contacts at z 0: 61 when
@ φ
′
@ t
′
=
1: 53 corresponding t o a g round surf ace settlement v elocity of
@ φ
′
@ t
′
Δ φ Δ z= Δ t =
0: 08cm=s . Liquef action occurs at t his location no t because of special prop-
erties of t he sand at t his point as all sand is treated equall y , but r at her be-
45
4.3. Qualitativ e R esults
0
0.5
1
1.5
2
2.5
0 0.2 0.4 0.6 0.8 1
Dimensionless Pressure
Dimensionless Height Above Datum
Total stress
Water Pressure
0
0.5
1
1.5
2
0 0.2 0.4 0.6 0.8 1
Dimensionless Perm.
Dimensionless Height Above Datum
F igure 4.1: P ar abolic w ater pressure i s tang ent t o t he linear t o tal s tress cur v e
at t he point of initial liquef action in a unif or m sand
cause of t he flo w and its inter action wit h t he boundar y conditions as will be
explained later . F or a 30 second eart hquak e and a 10 m deposit, t his im plies
a v ertical settlement of about 2.4 cm. N or mall y if t he g round drops se v er al
centimeters w e w ould see serious consequences f or pipelines or s tr uctures
on spread f oo tings.
Exam ple 2: Fine-scale v ar iability in per meability around an a v erag e v alue
Due t o t he im portance of per meability v ariation, w e consider a model of a
loose h ydr aulic fill sand whose la y ering s tr ucture causes a fine-scale oscilla-
tion in per meability of 30% about a typical v alue k
′
d
= 1+ 3= 10sin( 2π 40z
′
) .
46
4.3. Qualitativ e R esults
0
0.5
1
1.5
2
2.5
0 0.2 0.4 0.6 0.8 1
Dimensionless Pressure
Dimensionless Height Above Datum
Total stress
Water Pressure
0
0.5
1
1.5
2
0 0.2 0.4 0.6 0.8 1
Dimensionless Perm.
Dimensionless Height Above Datum
F igure 4.2: Liquef action initiates at essentiall y t he same place as in t he con-
s tant per m eability exam ple despite a fine-scale oscillation in per meability of
30%. The dominant ter m in t he differential equation is t he k
′
d
@
2
P
′
@ z
′ 2
ter m so
long as k
′
d
̸ 0 im pl ying t hat lar g e s table g r adients in pressure do no t f or m
unless per meability becomes quite small as seen in t he silt-la y er exam ple.
This oscillation occurs on a small lengt h scale O( 1= 40) of t he t o tal lengt h, and
is a sim ple model f or v ariability in h ydr aulic fill material and t he natur al seg-
reg ation of sandy soil b y v ar ying sinking r ates due t o v ar ying g r ain sizes as
pointed out b y (K okusho and K ojima 2002 ). N o te t hat t he w ater pressure
cur v e f or t his case is v er y similar t o t hat of figure 4.1 .
This exam ple allo w s us t o highlight t he role of t he v ariability in per me-
ability . If t he per meability field is O( 1) t hroughout, but has a fine scale os-
47
4.3. Qualitativ e R esults
cillation of dimensionless lengt h scale λ ≪ 1 t hen t he ter m
@ k
′
d
@ z
′
becomes
O( 1= λ) which can be quite lar g e. T o balance t his w e can eit her g ener ate cur -
v ature of t he pressure field t o induce diffusion, or if diffusion is quenched
b y k
d
@
2
P
′
@ z
′ 2
0 t hen(
@ P
′
@ z
′
+ 1) mus t be O( λ) t o k eep
@ k
′
d
@ z
′
(
@ P
′
@ z
′
+ 1) = O( 1) .
Suppose t hat a fine scale v ariation in P
′
de v elops so t hat P
′
= F( z
′
= λ)+
G( z
′
) where F is an O( 1) function representing t he fine scale v ariation and G
is O( 1) representing t he lar g e scale v ariation. Thenj
@
2
P
′
@ z
′ 2
j= O( 1= λ
2
) wit h a
f act or 1= λ coming from each deriv ativ e of F wit h respect t o z
′
. In or der f or
diffusion t o be quenched w e mus t ha v e k
′
d
≪ λ
2
so t hat k
′
d
@
2
P
′
@ z
′ 2
0 in equation
4.19 . Since in real materials 0 < k
′
d
< O( 1) , necessaril y t he region where
bo t h k
′
d
≪ λ
2
andj
@ k
′
d
@ z
′
j= O( 1= λ) is extremel y small. T o see t his, consider
t hat if k
′
d
is small, and t he deriv ativ e of per meability is neg ativ e, t hen it mus t
become nonneg ativ e quickl y t o k eep k
′
d
> 0 ; or if t he deriv ativ e is positiv e,
k
′
d
increases r apidl y so it does no t remain small f or long.
If t he o v er all size of k
′
d
is small e v er ywhere, t hen w e can rescale t he k
′
d
v ariable e v er ywhere, and b y rescaling t he t
′
v ariable t he equation remains
symmetric t hough no w all processes tak e long er in dimensional time. This
sho w s t hat f or lo w er per meability soils an o v er all lo w er le v el of f orcing is
required f or liquef action as t he same chang e in φ
′
w ould liquefy t he soil if
im posed o v er a long er dimensional time. In g ener al, t he conclusion is t hat
lar g e s table g r adients in pressure occur primaril y in small regions of v er y
lo w per meability soils surrounded b y higher per meability as sho wn in a later
exam ple wit h a t hin silty la y er .
By using a F ourier series approach t o appro ximate t he anal ytical solution
w e find t hat t he oscillations in per meability are com pletel y dominated b y t he
tendency f or pressure diffusion in t he s teady s tate case so t hat t he pressure
field is near l y indis tinguishable from t he ear lier cons tant per meability field
figure 4.1 , and does no t de v elop an im portant fine scale oscillation since t he
per meability ne v er become s small enough; in t his exam ple
@ φ
′
@ t
′
= 1: 45 and
z
′
= 0: 599 com pared t o t he cons tant per meability model where
@ φ
′
@ t
′
= 1: 53
and z 0: 61 . Since t he figure f or t his exam ple is v er y similar t o figure 4.1 , it
is reproduced in t he supplement.
Exam ple 3: Silty Inter -La y er Ano t her case of interes t is a sim ple model f or
t he experiment of K okusho ( 1999 ). In t his situation w e ha v e loose sandy soil
wit h k
′
d
= 1 and a t hin la y er of much lo w er per meability k
′
d
: 001 wit h
a lar g e g r adient in per meability around t his minimum (see figure 4.3 ). Be-
48
4.3. Qualitativ e R esults
cause of t he lar g e per meability g r adients, t here can be f as t changing bound-
ar y la y ers in which t he pressure g r adients are high and de v elop o v er small
regions so t hat t he cur v ature of pressure is quite lar g e. If w e sol v e t he s teady
s tate equation w e find t hat wit h
@ φ
′
@ t
′
= 1: 23 t he confining of pressure b y t he
inter -la y er f orces a localized loss of g r ain contacts at z
′
= 0: 497 jus t belo w
t he center of t he la y er . In t he experiment of K okusho ( 1999 ) t here is t he f or -
mation of a t hin w ater film belo w t he silt inter -la y er in a manner consis tent
wit h t hese results.
A f ollo w up t o t he abo v e experiment considered a silt inter -la y er in a fine
sand, in which a long las ting sable la y er of w ater belo w t he silt la y er w as
obser v ed (K okusho and K ojima 2002 ). The w ater la y er will be long las ting if
t he upper soil can no t f all quickl y t hrough t he w ater due t o dr ag. This can be
explained in g ener al b y t he r atio of dr ag f orce t o w eight f or small particles.
If w e consider a single spherical soil particle of r adius R , t he w ell kno wn f or -
mula f or S t ok es dr ag is f
d
= 6π μ R v (Brinkman 1949 ) whereas t he w eight
f or a spherical g r ain is w =
4
3
π R
3
ρ g . If w e tak e t he fluid v elocity t o be t he v e-
locity in t he w ater la y er immediatel y abo v e t he sand v=
k
d
μ
@ P
@ z
t he char acter -
is tic v alue of f
d
= w f or one nondimensional unit of each v ariable is t heref ore:
9k
d0
2G
g r
R
2
. F or sand g r ains wit h R 0: 5mm t his r atio is 0: 001 whereas f or a silt
g r ain wit h R 0: 01mm t he r atio is 2: 5 indicating t hat silt particles which
f all do wn int o t he w ater la y er are easil y pushed back up b y t he dr ag caused
b y fluid flo w whereas t he sand will easil y sink under t hese conditions. In
addition t o t his basic inter pretation, Brinkman ( 1949 ) and V erschoor ( 1951 )
sho w t hat t he dr ag on t he particle increases dr amaticall y as it integ r a tes int o
t he porous medium at t he w ater -silt boundar y so t hat once a la y er is f or med,
t he particles abo v e are unlik el y t o detach and sink. The effect of g r ain sizes
is discussed more com prehensiv el y in an additional section belo w .
In ano t her experiment b y K okusho and K ojima ( 2002 ), which in v ol v ed
a fine sand o v er l ying a coarse sand, a tr ansit or y turbulent w ater la y er w as
obser v ed at t he interf ace betw een t he sands. The differences betw een t his
coarse and fine sand case com pared t o t he fine sand and silt case can be
unders t ood in ter ms of tw o f act ors: t he lack of a long-ter m s trong g r adient in
pressure which w ould induce a lar g e flo w v elocity , and t he lar g er g r ain sizes
in v ol v ed in t he coars and fine sand experiment. In furt her cases, in v ol ving a
fine sand inter -la y er in a coarse sand, t he results w ere subs tantiall y similar t o
t he sand wit h silt inter -la y er experiment wit h an altered timescale due t o t he
different per meabilities in v ol v ed (K okusho and K ojima 2002 ). F inall y , t he
49
4.3. Qualitativ e R esults
aut hors alter t heir fine sand o v er l ying a coarse sand experiment where onl y
tr ansient turbulent w ater films f or med at t he interf ace b y reducing t he height
of t he w ater table. In t his case t he w ater la y er at t he interf ace is s tabilized
(K okusho and K ojima 2002 ). The lack of w ater abo v e t he w ater table means
t hat t he effectiv e size of t he upper fine sand la y er is onl y 7 cm com pared t o 90
cm in t he original experiment, since t he unsatur ated sand abo v e t he w ater
le v el does no t affect t he flo w and t he fluid pressure drops t o atmospheric
at t he fluid/air boundar y . This effectiv el y t hinner upper soil la y er means a
lar g er pressure g r adient can f or m and higher fluid v elocities de v elop, t hus
higher dr ag de v elops on t he upper sand particles which s tabilizes t he w ater
film.
Exam ple 4: Liq uefiable Sand Inter -La y er A typical w ar ning sign f or Geo tech-
nical engineers is t he exis tence of a loose sand la y er betw een l a y ers of denser
and especiall y lo w er per meability material. This is t he in v erse of t he silty
inter -la y er exam ple giv en abo v e. When w e mak e per meability k
′
d
= 1= 8 out-
side of a region of widt h about 20cm (10 times as wide as t he silty inter -la y er)
where per meability is 1 w e find t hat t he contr action of t he deposit pushes
w ater upw ar ds out of t he coarse sandy la y er and t he deposit liquefies in t he
fine sand abo v e at z
′
= 0: 65 . Also
@ φ
′
@ t
′
= 0: 198 which is muc h lo w er loading
t han in eit her t he unif or m or silty la y er case. In t his case, t he location and
manner of liquef action is less im portant t han t he magnitude of t he loading
required (see figure 4.4 )). Due t o t he lo w er per meability of t he main de-
posit, t he diffusion is less dominant, and wit h less dr ainag e, it is easier t o
cause liquef action. Ho w e v er , since t he g r adients in pressure remain O( 1)
t he perturbation ter ms do no t become im portant here. This highlights t he
phenomenon t hat “dr ainag e” is a non-local phenomenon which is affected
b y t he entire per meability field. Sand itself ma y be locall y w ell dr ained, but
can be confined b y lo w per meability boundaries.
Exam ple 5: Sand o v er lain b y silt This model (figure 4.5 ) is modeled af ter
t he centrifug e results of F ieg el and K utter ( 1994 ) in which a silt o v er lies a
loose sand, and is subject t o shaking in a g eo technical centrifug e. In t hese
experiments, pore w ater pressure w as recor ded in se v er al locations wit hin
t he sand and silt. F or our model of t his case w e ha v e rescaled t he z axis t o
correspond t o t he dep t h of t he experiment and correspondingl y rescaled P
′
,
φ
′
and t
′
. Contr activ e f orcing
@ φ
′
@ t
′
is spatiall y v ar ying and appreciable onl y in
50
4.3. Qualitativ e R esults
0
0.5
1
1.5
2
2.5
0 0.2 0.4 0.6 0.8 1
Dimensionless Pressure
Dimensionless Height Above Datum
Total stress
Water Pressure
0
0.5
1
1.5
2
0 0.2 0.4 0.6 0.8 1
Dimensionless Perm.
Dimensionless Height Above Datum
F igure 4.3: Silt la y er embedded wit hin loose sand. Liquef action in itiates be-
lo w t he region of lo w es t per meability where t he lar g e g r adient in per me-
ability induces a lar g e second deriv ativ e of pressure, and a t hin w ater film
f or ms.
51
4.3. Qualitativ e R esults
0
0.5
1
1.5
2
2.5
0 0.2 0.4 0.6 0.8 1
Dimensionless Pressure
Dimensionless Height Above Datum
Total stress
Water Pressure
0
0.5
1
1.5
2
0 0.2 0.4 0.6 0.8 1
Dimensionless Perm.
Dimensionless Height Above Datum
F igure 4.4: A high per meability la y er acts as a dr ain f or w ater pressure in
t he lo w er region and a source f or w ater pressure in t he upper region. Lique-
f action occurs abo v e t he high per meability la y er . This result occurs whet her
w e use unif or m f orcing t hroughout or localized f orcing onl y wit hin t he high
per meability la y er .
52
4.3. Qualitativ e R esults
t he lo w er sand region, t he silt region is assumed no t t o contr act, and hence
t he w ater pressure cur v e is near l y linear in t he silt which im plies fluid flo w
at cons tant v elocity . Clear l y t hen, w ater is squeezed out of t he lo w er sand
and pressure diffuses int o t he upper silt, liquefying it at t he base. The ex-
perimental results of F ieg el and K utter ( 1994 ) qualitativ el y ag ree wit h t his
s teady s tate solution. N o te t hat t he consequences of liquef action, a f ailure
via a sand-boil, w ere obser v ed t o occur much later , presumabl y near a lo-
calized w eakness wit hin t he silt la y er . In our exam ple, liquef action occurs
at z
′
= 0: 695 which is slightl y abo v e t he center of t he per meability tr ansi-
tion, which is located at z = 0: 629 . The required loading is
@ φ
′
@ t
′
= 0: 011
which ultimatel y is unrescaled (bo t h φ
′
and t
′
are rescaled b y t he same f ac-
t or , hence canceling). This is significantl y smaller t han in t he unif or m sand
exam ple (figure 4.1 ) where
@ φ
′
@ t
′
= 1: 53 . The relativ el y t hick silty la y er mak es
dr ainag e difficult as expected, and t his contributes t o a lar g e reduction in t he
necessar y loading.
Exam ple 6: Theoretical process of t her mal liq uef action of silt By subs ti-
tuting k
′
d
= ε
^
k
d
in t he incom pressible equation 4.19 w e can consider t he case
when typical per meability is small relativ e t o sand, such as ε < : 001 . When
t his is t he case, diffusion does no t act on t he same time scale, and in f act
t he diffusion time scale Δ t= ε > 140s is lar g er t han t he typical dur ation of
t he eart hquak e. This rescaling process also mak es t he γ
vol
ter m more dom-
inant, and t he full com pressibility effects might become im portant t hough
w e ignore t hem in t his exam ple as our pur pose is sim pl y t o highlight a ne w
phenomenon which could pro vide an a v enue f or furt her s tudy . If a silt is
someho w precluded from sus tained g r ain sk elet on contr action
@ φ
′
@ t
′
, per haps
due t o its initial porosity being near t he equilibrium s tate, it could s till t he-
oreticall y liquefy b y absor p tion of w a v e ener gy int o t her mal ener gy . The
sim ples t case is a single homog eneous la y er of silt similar t o t he homog e-
neous la y er of sand abo v e. Ins tead of contr activ e f orcing, w e appl y
@ S
′
@ t
′
and
attem p t t o liquefy t he silt via heating. The equation is similar t o unif or m
contr action, and produces t he same pressure field as figure 4.1 . W e find
t hat t he homog enous la y er liquefies at z
′
0: 61 when
@ S
′
@ t
′
= 1: 53ε= γ
S
. The
ques tion arises t hen as t o whet her ε= γ
S
= ord( 1) in some regimes of shak -
ing and g r ain size im pl ying t hat silt might t her mall y liquef y b y absorbing
w a v e ener gy t o heat. This is an interes ting ques tion f or which a sim ple cen-
trifug e experiment using a w ell insulated buck et and measuring bef ore and
53
4.3. Qualitativ e R esults
0
2
4
6
8
10
0 0.2 0.4 0.6 0.8 1
Dimensionless Pressure
Dimensionless Height Above Datum
Total stress
Water Pressure
0
0.5
1
1.5
2
0 0.2 0.4 0.6 0.8 1
Dimensionless Perm.
Dimensionless Height Above Datum
F igure 4.5: A model f or t he experimental results of F ieg el and K utter ( 1994 ) in
which t he s tep in per meability induces liquef action at t he sand-silt boundar y
where a shar p kink in pressure occurs. Here t he spatial dimension has been
rescaled b y a f act or ε : 34 so t hat z
′
2 [ 0; 1] corresponds t o z 2 [ 0; 3: 5]m ,
t he dimensions of t he experiment. Pressure, φ
′
, and time ha v e been rescaled
b y ε
2
0: 12 t o maintain t he or der of magnitude of t he coefficients in t he
equation. Scaling of t he f orcing
@ φ
′
@ t
′
cancels out s o t hat numerical v alues are
directl y com par able t o o t her exam ples. F orcing is cons tant in t he sand, and
drops t o near zero jus t belo w t he per meability tr ansition, hence onl y t he sand
contr acts appreciabl y . The g r adients s ta y small enough t hat t he perturbation
ter ms do no t become im portant.
54
4.4. The role of g r ain sizes and dr ag f orces f or liquef action consequences
af ter tem per atures could giv e insight, or per haps borehole measurements of
t her mal chang es during an eart hquak e might be a v ailable t o pro vide insight.
4.4 The role of g rain sizes and drag f orces f or
liq uef action conseq uences
Since w e are considering t he liquef action of silt, it should be mentioned t hat
pure silts and pure g r a v els are g ener all y t hought t o be lo w er risk f or liquef ac-
tion t han sands, t hough a more nuanced vie w of silt and g r a v el liquef action
is de v eloping (cf. (Idriss and Boulang er 2006 ; Seed et al. 2003 )). Examining
t he cases f or different g r ain sizes wit hin t he context of sim ple s teady s tate
w ater flo w described here can be ins tr uctiv e.
In or der f or sus tained contr action of t he g r ain sk elet on t o occur w e need
g r a vitational and o v er -bur den f orces on t he particles t o be balanced b y dr ag
f orces of t he w ater so t hat t he g r ains tr a v el do wn w ar ds at a slo w ter minal
v elocity while t he w ater is pressurized and tr a v els upw ar ds t o w ar ds t he sur -
f ace. Because of t his balance of f orces, and t he r apidity of pressure diffusion,
t he details of g r ain mo tion at each point in t he loading cy cle are less im por -
tant t han t he o v er all a v er ag e contr activ e trend.
In t he exam ples b y K okusho and K ojima ( 2002 ), t he sand particles de-
scend r apidl y t hrough t he w ater while t he silt particles sit on t op of t he w ater
la y er . These settlement v elocities are related t o t he dr ag f orces on t he parti-
cles. Ho w e v er when all particles are silt sized, e v en if all particle contacts are
los t, settlement will be extremel y slo w and particles will no t mo v e f ar from
t heir initial positions during a loading cy cle due t o high dr ag f orces relativ e
t o g r a vitational settlement f orces. Consequences of such an e v ent on le v el
g round ma y be significantl y less des tr uctiv e t han in sand where dr ag f orces
do no t affect mobility as s trongl y and settlement is more r apid. F or g r a v els
mobility is e v en more r apid t han sand, but contr action of t he g r a v el does no t
induce s trong pressurization since t he pressure dissipates extremel y r apidl y
due t o t he high per meability , unless surrounding lo w er per meability soils
im pede flo w . This im plies w e should expect f as t settlements but onl y tr an-
sient liquef action in contr activ e g r a v els.
Alt hough in t he v ertical direction w e expect r apid settlement from g r a v el
particles and slo w settlement from silt, in ter ms of shear displacements t he
roles ma y be re v ersed. The r apidity wit h which g r a v el particles descend
55
4.5. The problems wit h small scale em pirical tes ts
ma y allo w t hem t o reg ain contact more r apidl y t hereb y reducing t he time
during which shear s trengt h is los t, whereas silt particles will tak e a long
time t o reg ain t heir shear -tr ansf erring ability t hrough settlement. Soils wit h
a wide r ang e of particle sizes will be more com plex. The de v elopment of a
soil beha vior model which tak es t he fluid dr ag f orces int o account explicitl y
is a w ort h y goal. Whate v er t he case, t his anal ysis highlights t he g r adient in
particle size, which is also related t o t he g r adient in per meability , as pla ying
an im portant role in t he liquef action process especiall y wit h reg ar ds t o t he
de v elopment of fluid lubr ication la y ers as seen in K okusho ( 1999 ); K okusho
and K ojima ( 2002 ).
4.5 The problems wit h small scale em pir ical tests
T o unders tand t he role of w ater flo w on small scale sam ples, consider t he ex-
am ple in figure 4.1 where all t he sand is unif or ml y similar in initial porosity
and per meability conditions. If w e tak e a slice of dimensionless lengt h 0.01
from t his exam ple, it corresponds t o a sam ple similar in size t o a h ypo t heti-
cal triaxial sam ple of about 10cm char acteris tic lengt h. R ecall t hat all points
in t he model of exam ple 1 w ere squeezed equall y since
@ φ
′
@ t
′
w as a cons tant
in space. The r ate of squeezing in exam ple 1 w as 1.53 units of Δ φ per 1 unit
of dimensionless time, or a v olumetric s tr ain of about 1% in 130 s which is
per haps a reasonable timescale f or a tablet op undr ained liquef action exper -
iment. In t he absence of dr ainag e, at t his v olumetric s tr ain r ate, w e w ould
g ener ate in our h ypo t hetical sam ple 100kP a of additional pressure e v er y 0.6
seconds, and at t he end of 130s w e w ould ha v e about 22.5 MP a of w ater
pressure. Ins tead, in exam ple 1, w e ha v e about 100 kP a of ex cess pressure
(1 dimensionless unit o v er atmospheric) at infinite time, as sho wn in figure
4.1 at z = 0: 6 where liquef action occurs. This discrepancy is due t o dr ainag e
in s teady s tate conditions. N o te t hat t he pressure field in exam ple 1 is in-
fluenced b y t hree f act ors: t he r ate of squeezing (
@ φ
′
@ t
′
= 1: 53 ), t he cons tant
per meability , and t he boundar y conditions (im per meable at z = 0 , cons tant
atmospheric pressure at z = 1 and zero later al flo w). This is in contr as t t o t he
triaxial tes t where dr ainag e is pre v ented (im per meable at t op and bo tt om),
and t he boundar y conditions are much closer , and more com plex, due t o t he
surrounding r ubber membr ane and cons tant pressure bat h.
In exam ple 1, t he reason t he point at z = 0: 6 has liquefied is no t t hat it
has been squeezed har der , or t hat its g r ain sk elet on w as more suscep tible
56
4.5. The problems wit h small scale em pirical tes ts
t han o t her points, it is t hat fluid is flo wing continuousl y t hrough t he sam-
ple, and due t o a combination of squeezing, and imbalance in fluid flo w , t he
density of w ater at t his point is jus t slightl y higher t han it w as bef ore w e be-
g an squeezing t he sam ple. The role of w ater flo w is highl y nontrivial e v en
in t his extremel y sim ple exam ple.
In contr as t, consider taking our slice of dimensionless lengt h 0.01, re-
mo ving it from t he g round, and placing it int o a triaxial machine. N o w ,
t he boundar y conditions ha v e been brought 100 times closer t o t he sam ple
(about 10cm i ns tead of 10m) and ha v e been chang ed. In a triaxial tes t, w e
ha v e an im per meable but s tretch y membr ane surrounded b y a cons tant pres-
sure on t he sides, wit h a rigid base and loading cell belo w and abo v e. This
can be com pared t o t he la y ered soil in our mat hematical model where lat-
er al symmetr y im plies zero flo w horizontall y , and w e ha v e an im per meable
la y er belo w and cons tant atmospheric pressure abo v e. Alt hough a triaxial
tes t is no t reall y a 1D sys tem, w e can begin t o qualitativ el y anal yze t he ef-
f ect of in v es tig ating t he small scale sam ples using t he same equation 4.19 b y
choosing ne w scale cons tants t o inter pret t he results.
W e lea v e t he char acteris tic pressure at P
0
= 100KPa since t his is t he or -
der of magnitude of pressure w e mus t induce t o liquefy t he sam ple. Since
Δ φ depends onl y on t he pressure scale, it remains unchang ed. W e ignore
heating f or t hese pur poses. W e are lef t wit h re-e v aluating Δ z and Δ t . The
char acteris tic size of our sam ple is Δ z = 0: 1m , and t he ne w v alue of Δ t is
deter mined from Δ z and equal t o Δ t =
( Δ z)
2
φ
0
μ
k
d0
K
B
1: 3 10
5
s or 13 microsec-
onds. The triaxial boundar y conditions are no w com plex because t he y in-
v ol v e t he inter action of a pressure bat h wit h a r ubber membr ane which t hen
inter acts bo t h wit h t he g r ains and t he pore w ater , and t he y are no long er 1D.
Ho w e v er , it is saf e t o sa y t hat t he pressure ex erted b y t he membr ane on t he
pore w ater is relativ el y unif or m. N ext w e consider t he loading,
@ φ
′
@ t
′
. It is saf e
t o sa y t hat t his is near l y zero on our ne w scales, as t he loading is such t hat
t he or der of magnitude of pressure chang es t hat t hese experiments induce
is enor mousl y less t han P
0
= 100kPa in 13 microseconds.
In t he absence of an y significant loading, t he w ater pressure inside t he
sam ple is in g r a vitational equilibrium at all times, since t he equilibriation
time is 10 million times shorter t han t he dur ation of t he experiment, and
since our ne w lengt h scale is Δ z = 10cm t he g r adients in pressure are tin y
as w e will ne v er induce an ywhere near 100 kP a fluid pressure differences
across t he 10cm sam ple.
57
4.6. Conclusions
Jus t as in t he full scale equation, t he w ater pressure in our im per meable-
membr ane-wr apped triaxial sam ple is dependent on t hree t hings: t he dy -
namics of t he loading which w e ha v e sho wn t o be essentiall y zero, t he uni-
f or m per meability which remains unchang ed, and t he pressure induced b y
t he boundar y conditions which in v ol v e bo t h g r ains and w ater inter acting
wit h a highl y com pliant im per meable m embr ane surrounded b y a cons tant
pressure bat h. In t he triaxial tes t t he fluid pressure is entirel y deter mined
b y t his boundar y condition as all points are in equilibrium wit h t he bound-
ar y . T o t he extent t hat w e are no t interes ted in ho w r ubber membr anes inter -
act wit h g r ains and w ater under cons tant exter nal pressure conditions, w e
should also be highl y sk ep tical of an y results from tablet op triaxial or similar
specimens in so f ar as t he y are supposed t o represent ho w t he soil beha v es
in t he g round. As w e ha v e sho wn at t he beginning of t his section, under
t he g round, w ater flo w is extremel y im portant, and tablet op experiments,
as currentl y perf or med, sim pl y can no t model t his process.
In contr as t t o triaxial type experiments, g eo technical centrifug e experi-
ments such as F ieg el and K utter ( 1994 ) ins tead use scaled fluid viscosity so
t hat t he diffusion time from a point deep in t he model t o a point on t he
surf ace where cons tant atmospheric pressure exis ts will scale appropriatel y .
Also t he model is built in la y ers, so t hat t he per meability field of t he entire
deposit is modeled correctl y , including per meability g r adients. Theref ore
fluid flo w can be proper l y accounted f or in g eo technical centrifug es t hrough
t he proper application of scaling la w s.
4.6 Conclusions
By deriving a com prehensiv e equation f or t he w ater pressure de v elopment
wit hin la y ered soils and anal yzing some critical exam ples w e ha v e sho wn
t hat fluid flo w driv en b y porosity chang es and per meability chang es is t he
main f act or t hat causes pore pressure chang e during liquef action. By com-
parison wit h experimental results (F ieg el and K utter 1994 ; K okusho 1999 ;
K okusho and K ojima 2002 ) w e sho w t hat sandy soil is w ell dr ained enough
t hat t his soil comes t o s teady s tate on t he same time scale as t he eart hquak e.
Also, due t o our model of t his dr ainag e, w e can identify se v er al modes of liq-
uef action, from fluid flo w induced “f as t settlement” t o localized w ater film
f or mation at per meability boundaries, t o flo w induced liquef action of lo w
per meability la y ers t hat o v er l y a high per meability la y er , and e v en po ten-
58
4.6. Conclusions
tiall y t her mall y induced liquef action when per meability is v er y lo w . In each
of t hese exam ple cases, w e do no t attem p t t o couple t he beha vior of t he g r ains
t o t he beha vior of t he w ater , or pro vide time-v ar ying dynamic solutions, as
no lar g e scale continuum model is a v ailable f or g r ain-fluid inter action be-
ha vior com par able in quality t o t he DEM model of Goren et al. ( 2011 ). Our
intention in t his paper has been t o highlight t he im portance of fluid flo w and
spur an unders tanding of t he need t o de v elop soil models t hat can account
f or flo w .
Modeling a 1/4 m square cross section 10 meters long wit h 1mm diame-
ter typical g r ain sizes w ould require per haps 10
8
g r ains and be com putation-
all y prohibitiv e. Ho w e v er f or qualitativ e unders tanding of t he liquef action
process it is sufficient t o perf or m a careful anal ytical anal ysis of a contin-
uum model of w ater flo w t hrough soil driv en b y chang es in pore v olume.
Our anal ysis pro vides tremendous insight wit hout supercom puting require-
ments, from which w e mak e se v er al im portant conclusions:
• Soil liquef action in w ater satur ated sands is a process usuall y in v ol v -
ing fluid flo w due t o g r a vitational pum ping, no t purel y fluid s tagna-
tion in an “undr ained” e v ent. Thin regions of lo w per meability silt do
no t pre v ent dr ainag e, as lar g e pressure g r adients will de v elop t o f orce
fluid t hrough t hese la y ers. Ho w e v er , t heir presence significantl y in-
fluences t he location and ease of onset of liquef action. Lar g e regions
of lo w per meability soil can significantl y reduce dr ainag e, leading t o
liquef action at significantl y lo w er r ates of g r ain contr action.
• Due t o w ater pressure diffusion, w e can no t examine t he liquef action
process b y sim pl y considering local properties of soil in a small sam-
ple on a tablet op appar atus, giv en t he char acteris tic lengt h scales of
ord( 10) meters, and corresponding char acteris tic time scales of tent hs
of seconds. In a tablet op appar atus wit h sam ples on t he or der of 10cm
t he timescale f or diffusion of pressure t o t he im per meable r ubber bound-
ar y and hence t o t he cons tant pressure boundar y conditions is on t he
or der of 10 microseconds, so t hat t he pressure t hroughout t he sam-
ple is equal t o t he boundar y pressure at t he membr ane. Theref ore,
g eo technical centrifug e experiments are t he onl y meaning ful ph ysical
experiments f or t his phenomenon.
• In lo w er per meability materials such as t hick regions of fine sand or
silt, t he time scales f or t he w ater flo w process can be long er . These long
59
4.6. Conclusions
char acteris tic times occur when t he typical per meability t hroughout
t he material isO( 10
13
)m
2
which is w ell belo w t he sandy r ang e (Bear d
and W e y l 1973 ). Due t o g r ain size effects, fluid flo w in silt causes high
dr ag on t he particles, limiting t heir mobility . Liquef action in extensiv e
silt deposits is an interes ting and different flo w regime.
• Im pedance t o flo w caused b y local shar p reductions in per meability is
t he main driving f orce in localized liquef action e v ents, and can cause
t he f or mation of w ater films t hat lubricate la y ers of soil allo wing lar g e
later al spreading. F or le v el g round, f ailure is expected via localized
sand boils at fissures, or t he liquef action of o v er l ying la y ers as fluid
flo w s tagnates. This phenomenon can onl y be explained b y explicitl y
ackno w ledging t he flo w of w ater .
• “Liquef action” as t he ter m is commonl y used is no t necessaril y asso-
ciated wit h t o tal loss of effectiv e s tress as is already kno wn. Lar g e set-
tlements could occur wit hout t he f or mation of a full y liquefied la y er if
loose soil settles while w ater flo w s t o t he surf ace. This ma y especiall y
occur in contr activ e g r a v el or coarse sand fills.
• In silts, where per meability is small and diffusion can no t act on t he
eart hquak e dur ation timescale, extensiv e shaking might concie v abl y
cause liquef action t hrough heating effects, a pre viousl y unexplored
phenomenon.
60
Chapter 5
Dissipation of W av es in a M olecular T est
Problem
What ar e t he wild w a v es saying
Sis t er , t he whole day long,
That ev er amid our playing,
I hear but t heir low lone song?
— Joseph Edw ar ds Car penter
5.1 Introduction
When t hinking about ho w t o model t he inter actions betw een g r ains and
eart hquak e w a v es which lead t o t he dissipation o f w a v e ener gy and rear -
r ang ement of g r ains it became clear t hat s tepping directl y t o t his com pli-
cated problem w as prohibitiv e. Ins tead, during discussions wit h eart hquak e
ph ysicis ts a sim pler tes t problem emer g ed. In t his problem w e attem p t t o
model t he cause of dissipation in w a v es tr a v eling t hrough a molecular pol y -
cr ys talline medium, which can be simulated t o high precision using MD
techniques. In t he absence o f w a v e spreading, exter nal r adiation, and scat-
tering, where all ener gy remains wit hin t he bar , an y w a v e dissipation mus t
be described via some o t her mechanism. The results of t his in v es tig ation are
reproduced belo w and are being prepared in subs tantiall y similar f or m f or
furt her publication.
61
5.2. The Problem
5.2 The Problem
A t hin rod wit h a square cross section made of pol y cr ys talline molecular ma-
terial sits in free space and t he tw o ends are tapped t o induce tw o opposite-
tr a v eling com pression w a v es. The net momentum and angular momentum
of t he rod is zero and w e assume no mechanism f or t he t o tal ener gy of t he
sys tem t o chang e, such as acous tic or t her mal r adiation or o t her coupling t o
t he surroundings. The classic w a v e equation deriv ed from linear elas ticity
is:
@
2
u
@ t
2
=
@ v
@ t
= v
2
w
@
2
u
@ x
2
(5.1)
Where v=
@ u
@ t
is a v elocity field along t he lengt h of t he bar , v
w
=
√
E
ρ
is
t he w a v e v elocity and u is t he displacement field of t he material from equi-
librium position.
This classic equation, deriv ed from t he assum p tion of s tress proportional
t o s tr ain
@ u
@ x
predicts no dissipation whatsoe v er and t he bulk scale w a v es con-
tinue t o tr a v el back and f ort h f ore v er in such a model. Ho w e v er e v en wit hout
ha ving perf or med t he experiment, w e are f air l y sure t hat t here is some long-
time limiting beha vior in which t he rod ceases t o ha v e a discer nible bulk -
scale w a v e, and ins tead is sim pl y at higher tem per ature t han bef ore. The
w a v e equation is t heref ore t he asym p t o tic limit where obser v ational size is
much lar g er t han molecular size so t hat measurements are no t made of indi-
vidual m olecules, and dissipativ e effects are negligible on t he obser v ational
time scale. Depending on t he size of t he rod, t he tem per ature, and t he ma-
terial it is made of, t he time scale f or dissipation ma y be an yt hing from fr ac-
tions of a second t o per haps y ears. The ques tion is, ho w can w e use our
kno w ledg e of t he molecular nature of t he material t o produce predictions
f or t he bulk scale deca y of w a v es in a medium where dissipation is non-
negligible, and what are t he scaling la w s t hat deter mine when dissipation is
im portant?
In realis tic situations, in addition t o t he fundamental s of local entrop y
g ener ation, g eometric spreading, diffr action, and scattering of w a v es b y in-
homog eneities in t he material oper ate simultaneousl y . A t lar g e lengt h and
ener gy scales in a 3D medium, g eometric spreading o v er t he spherical w a v e
front, diffr action, and scattering com pletel y dominate o v er local entrop y g en-
er ation. When adding f or mal dam ping ter ms t o t he w a v e equation, t he choice
62
5.3. The Approach
of dam ping coefficient is usuall y made b y ref erence t o s tatis tical fitting of
obser v ational data in which t hese lar g er scale processes dominate. Such s ta-
tis tical fitting does no t giv e us a satisf act or y descrip tion of t he inherent local
entrop y g ener ation r ate. Onl y a model built on our kno w ledg e of t he micro-
scale inter actions can inf or m t he causal ques tions of ho w dissipation occurs
in t he absence of spreading and scattering.
5.3 The Approach
W e simulate a solid built of unif or m Lennar d-Jones (LJ) particles, inter acting
wit h a tr uncated po tential, using t he LAMMPS molecular dynamics simu-
lat or (Plim p t on 1995 ). U nits are t he non-dimensional “Lennar d-Jones” units
described in t he LAMMPS manual. In particular 3T= 2 is t he a v er ag e kinetic
ener gy per molecule at t her mal equilibrium, all molecules ha v e mass 1, and
1 unit of ener gy is t he escape ener gy required t o separ ate tw o molecules in
free space initiall y at equilibrium dis tance. Lengt hs are measured as mul-
tiples of t he zero-crossing lengt h f or t he po tential, and time is defined b y
m v
2
= 2 = 1 unit of ener gy when one molecule tr a v els at v elocity v equal t o
one unit of dis tance in one unit of time.
W e repeatedl y choose a r andom number of particles N , an initial tem-
per ature T
0
and a v elocity pulse s trengt h v
0
t o r un 120 different simulations.
The par ameters w ere chosen so t hat N = 2500 + n where n w as exponen-
tiall y dis tributed wit h a v er ag e 50000. Using LJ w alls, w e confine t he par -
ticles int o a region whose dimensions are appro ximatel y ( δ; δ; 10δ) in t he
x; y; z directions where δ = ( N= 10)
1= 3
. Af ter initiall y equilibriating a bar
at higher tem per ature, and remo ving t he boundar y w alls, w e slo w l y lo w er
t he tem per ature and bring t he bar t o a ne w s tate of equilibrium at tem per -
ature T
0
. W e choose T
0
2 [ 0: 025; 0: 25] unif or ml y so t hat molecules are no t
g ener all y ho t enough t o boil off t he surf ace of t he bar in an y sam ple simu-
lation. Then w e select t he molecules in a region of lengt h 3δ= 2 at t he tw o
ends of t he z axis of t he bar and superim pose a v elocity v
0
on t op of t heir
t her mal v elocity dis tribution, oriented t o w ar d t he center of t he bar . The v e-
locity pulse s trengt h w as chosen unif or ml y at r andom in[ 0: 05; 0: 25] . F inall y
w e zero out t he linear and angular momentum exactl y so t hat t he bar does
no t mo v e in space and w atch t he time e v olution of t he sys tem. A t each slice
of widt h δ= 2 along t he z axis, w e output t he number of molecules, t he net ki-
netic ener gy , t he net po tential ener gy , and t he center of mass v elocity wit hin
63
5.4. The Model
each slice at each recor ded times tep. W e define t he t o tal w a v e kinetic en-
er gy as W =
∑
N
slice
i= 1
n
i
v
2
i
= 2 . Due t o fluctuations in t he initialization process
N
slice
is no t alw a ys 20. F or exam ple, occasionall y a 21s t bo x might ha v e a f e w
molecules in it, but w e can idealize t he sys tem as nominall y ha ving N/20
molecules in each of 20 slices of widt h e xactl y δ= 2 along t he lengt h and t his
idealization is good t o wit hin a f e w percent.
The t o tal kinetic w a v e ener gy deca ys wit h time accor ding t o a consis tent
patter n of time e v olution which v aries wit h N , T , and v
0
. If w e assume a
perf ect linear elas tic w a v e equation t hen w e can define a ne w time scale b y
t
w
= ( L= 2)= v = ( L= 2)=
√
K= ρ . W e find experimentall y t hat ρ 1: 05 in t hese
sys tems, and assume t hat K =
@
2
U
@ r
2
j
r
min
= 9( 2
( 8= 3)
) 57: 15 leading t o a w a v e
speed of about 7: 38 in molecular units. The time scale t
w
is t he molecular
time required f or our w a v es t o tr a v el from t he ends of t he idealized bar t o
t he center . In perf or ming s tatis tical inf erence later w e will fit t he t
w
time scale
more precisel y .
W it h precisel y calibr ated timescale, w e t heref ore expect t he t o tal w a v e
kinetic ener gy t o reach minimum at times t= t
w
=f 1; 3; 5;:::g . If w e choose
a bar scale lengt h of l
w
= δ= 2 t hen t he w a v e speed in our ne w bar scale units
is nominall y 10 , and a sys tem which is purel y linear elas tic w ould collapse
do wn t o a single solution ha ving a w a v e tr a v eling from t he end t o t he center
in one unit of time and bouncing back and f ort h wit hin t he bar infinitel y .
If dissipation occurs ho w e v er , t he dissipation need no t be independent
of N in t he ne w time and lengt h units as dissipation ma y depend on N in a
different w a y , as w ell as on T or o t her properties. Our next task is t o deter -
mine ho w t his bar beha v es including dissipation as a function of N and T .
W e mention f or com pleteness t hat w e ignore t he effect of surf ace ener gy and
surf ace tension here, and find t hat it does no t mak e a lar g e difference in our
results.
5.4 The Model
The dynamics of our bar are represented at t he bulk scale as 20 measure-
ments of different properties, but especiall y t he number of particles and
center of mass v elocity of each slice. Of course it is possible t o choose an-
o t her number of bo x es, but so long as t he number is fix ed across t he dif-
f erent sam ples and is lar g e enough in e v er y sam ple t o contain a nontrivial
s tatis tical sam ple of particles t hen it mak es no difference what t he number
64
5.4. The Model
Box 1 Box 2 Box 0
/2
dx
s
F igure 5.1: An idealization of t hree slices of t he bar from which s tatis tical
a v er ag es are calculated in t he molecular dynamics simulation
is. What does matter ho w e v er is t he widt h of t he slices relativ e t o t he inter -
action lengt h of t he LJ po tential as w ell as t o t he f eature size of t he initial
w a v e pulse. The non-dimensional size of t he slice δ= 2 = ( N= 10)
( 1= 3)
= 2 and
non-dimensional tem per ature T are t heref ore t he par ameters of interes t f or
our dissipation model.
The deriv ation of our model proceeds accor ding t o a tem plate described
in chap ter 3 . W e s tart b y taking int o account t he molecular processes wit hin
a slice. The bulk scale a v er ag e s tatis tics are f or med from a v er ag es o v er t he
molecules in our slice. W e choose a slice on t he interior of t he bar and call
it slice 1, so t hat t he slice t o t he lef t is slice 0 and t he slice t o t he right is
slice 2 (see figure 5.1 ). The slice is fix ed in space in an Eulerian coor dinate
sys tem. Assume t hat t he bulk scale measurements at time t = 0 are n
i
; v
i
; T
i
.
W e choose a short time dt which is small com pared t o t he time t hat t he w a v e
65
5.4. The Model
tak es t o tr a v el t he widt h of t he slice δ= 2 and treat it as if it w ere a nons tandar d
dt in t he IS T model of nons tandar d anal ysis (N elson 1977 ).
The expected momenta P at times t , and t+ dt in bo x 1 are:
P
1
( t) = n
1
v
1
(5.2)
P
1
( t+ dt) = n
11
( v
1
+ dv
1
)+ dn
01
v
01
dn
10
v
10
dn
12
v
12
+ dn
21
v
21
(5.3)
Here n
11
ref ers t o t he number of particles t hat s tart in slice one and end
in slice one. The small chang es dn
i j
ref er t o t he number of particles mo ving
from slice i t o slice j , and t he v elocities v
i j
ref er t o t he mean v elocities of t hose
particles which s tart in slice i and end in slice j .
Since counts such as n
1
are O( 100) t o O( 10000) in our simulations it is
no t reasonable t o immediatel y ignore t he ir discrete nature and treat t hem as
continuous wit hout furt her anal ysis.
The onl y tr ul y continuousl y divisible fundamental quantities in t he sys-
tem are time and le ngt h. Deriv ed quantities momentum and ener gy are
t heref ore also continuousl y divisible, but bo t h also depend on t he quantized
mass
1
.
W e consider t he chang es in momentum on time scales t hat are short com-
pared t o t he time it tak es a molecule t o go a significant fr action of t he typical
inter -molecular dis tance (t h e μ dt time scale) and on long er time scales t hat
are s till short com pared t o t he time it tak es t he center of mass of t he bo x t o
tr a v el half t he bo x lengt h at its initial v elocity (t he dt time scale).
There are tw o contributions t o t he chang e in momentum, t he firs t is t he
quantized tr ansport of momentum across t he bo x boundar y b y tr ansport of
mass. N amel y when a molecule enters or lea v es t he bo x its momentum be-
comes part of or ceases t o be a part of t he a v er ag e o v er t he bo x. These tr an-
sitions in v ol v e jum ps of size proportional t o 1= n where n is t he number of
molecules in t he bo x, and are a P oisson-lik e process wit h t he particles be-
ing discrete, ho w e v er when n is reasonabl y sized t he jum ps are no t lar g e
relativ e t o t he a v er ag e, and t heref ore 1= n= dt = O( 1) so t hat w e can treat
t he momentum function as if it w ere an s-continuous nons tandar d function.
This “con v ectiv e deriv ativ e” portion of t he momentum tr ansport occurs at
1
The use of t he ter m “quantized” ref ers t o t he f act t hat mass comes in discrete molecules,
but all calculations are based on classical ph ysics wit hout an y of t he w a v e-function or o t her
aspects of quantum mechanics. The main issue is t hat due t o t he discrete masses, s tatis tical
a v er ag es o v er a bo x will be discontinuous when molecules cross bo x boundaries.
66
5.4. The Model
v elocities related t o t he center of mass v elocity of each slice so t hat t he dt
time scale is t he one of interes t.
Secondl y , t here is momentum tr ansport at infinite speed due t o t he f orce
fields from t he LJ po tential. W e do no t consider relativis tic effects and retar -
dation of t he f orces, especiall y at t he tin y dis tances modeled b y t he LJ inter -
action. Alt hough t his tr ansport is at infinite speed, t he r ang e is s till on t he
or der of a f e w molecular lengt hs. F orce-based momentum flux es t hrough
t he boundaries of t he bo x are t heref ore dominated b y t he inter actions near
t he surf ace of t he bo x wit hin a dis tance of dx
s
1 molecular lengt hs. In-
betw een jum p e v ents caused b y molecules crossing t he boundar y , t his mo-
mentum tr ansport is a continuous function, but due t o t he f act t hat t here are
multiple molecules in t he bo x and each one is mo ving at a different speed
and at a different dis tance from t he boundar y w e expect detailed fluctua-
tions in t his f orce which occur on t he μ dt time scale. F rom t he perspectiv e
of modeling t he increments at t he dt time scale each increment is itself an
integ r al o v er t he finer time scale. N ons tandar d numbers suppl y us wit h a
model of t his type of scale separ ation. If t
is one unit of non-dimensional
bar time, t hen dt = t
= I
1
and μ dt = dt= I
2
wit h I
1
; I
2
being nons tandar d in-
teg ers, and I
1
= I
2
= 0 is a model of t hree separ ate time scales each smaller
timescale infinitesimal wit h respect t o t he next lar g er one. W e t hink of t he
t
time scale as t he time it tak es t he w a v e t o mo v e half t he bar lengt h, dt as
t he time it tak es t he w a v e t o mo v e a microscopic bar lengt h, and μ dt as on
t he or der of t he time it tak es f or t he smalles t oscillation in t he f orce along a
boundar y due t o tin y r earr ang ements of t he molecules. Alt hough in reality
neit her I
1
nor I
2
are tr ul y nons tandar d, t he concep t of a nons tandar d integ er
is a good model f or a “lar g e but no t exactl y specified” number which aligns
w ell wit h t he asym p t o tic ideas w e use t o build t his model.
Momentum T ranspor t b y Mass T ranspor t
Alt hough t he anal ysis of t he momentum tr ansport b y particles crossing t he
slice boundaries is interes ting w e find t hat f or our pur poses it can be ignored
in t hese simulations. If w e let t he slice mo v e wit h t he center of mass, in a
Lag r angian fr ame of ref erence, t hen t he s tandar d equation connecting t he
r ate of chang e of t he Eulerian momentum a v er ag e P t o t he net f orces is:
@ P
@ t
+ v
@ P
@ x
= F
net
67
5.4. The Model
Since f or typical experiments, v = O( 0: 1) and n > 100 whereas dx > 5 ,
t he contribution of t he con v ectiv e portion v
@ P
@ x
< 0: 1= 500 < : 0002 whereas
w e can es timate t he or der of magnitude of t he typical f orce as t he chang e in
momentum from t he peak of t he w a v e t o zero as t he w a v e passes our slice,
F
net
= O( Δ P= Δ t) = n v=( 3dx= v
w
) = ( 100)( 0: 1)=( 15= 7: 4) = : 5 . Con v ectiv e
tr ansport of momentum is t heref ore a small fr action of t he t o tal tr ansport
and will be considered as part of t he noise f or sim plicity .
Momentum T ranspor t b y F orces
Consider a time increment in which no molecules cross an y boundaries, so
t hat all t he momentum tr ansport is via f orces. In a time dt t he net f orce on
t he molecules in bo x 1 can be described as:
F
1net
=
(
( ρ
1
+ ρ
2
)
2
F
l j
( r
e
12
)+
( ρ
0
+ ρ
1
)
2
F
l j
( r
e
01
)
)
( A dx
s
+ ε
n
)
Here A is t he effectiv e cross sectional area of t he bar ρ
n
is t he number
density of t he bar in bo x n , F
l j
is t he Lennar d-Jones f orce function, and r
e
i j
is
t he effectiv e dis tance betw een t he bo x es which giv es t he per -molecule f orce
needed t o account proper l y f or t he actual t o tal f orce. The v ariable dx
s
is a
short lengt h of ord( 1) molecular lengt h such t hat a v olume A dx
s
contains
t he v as t ma jority of t he molecules in bo x 1 t hat are inter acting appreciabl y
wit h molecules in t he neighboring bo x on each side. Our ε
v
is t he error in t he
inter action v olume which w e tak e t o be small enough t hat it can be modeled
as infinitesimal relativ e t o A dx
s
.
If w e f ollo w all t he molecules t hrough t heir com plex mo tions on t he μ dt
time scale as is done in t he molecular dynamics simulation, and calculate
t he effectiv e r v alues w e w ould see a continuous tr a ject or y of r
e
. Ho w e v er
w e wish t o mak e our model at a long er time-scale, one more appropriate f or
t he mo tion of t he center of mass of t he molecules in t he bo x. This dt time
scale is a lar g e multiple of t he μ dt time scale. F or t heoretical pur poses, t he r
e
v alues at t he end of our dt time s tep can be modeled as some drif t, which is
closel y related t o t he center of mass v elocity , and some r andom com ponent
which is related t o t he detailed wiggles of t he continuous pat h of r
e
. W e tak e
t hen:
68
5.4. The Model
dr
01
= ( v
1
v
0
) dt+ σ( T)◦ dQ
01
(5.4)
dr
12
= ( v
2
v
1
) dt+ σ( T)◦ dQ
12
(5.5)
ρ
01
= ( ρ
0
+ ρ
1
)= 2 (5.6)
ρ
12
= ( ρ
1
+ ρ
2
)= 2 (5.7)
dP
1force
=
(
ρ
01
( F
l j
( r
e
01
)+ F
l j
( r
e
01
+ dr
01
))
2
+
ρ
12
( F
l j
( r
e
12
)+ F
l j
( r
e
12
+ dr
12
))
2
)
A dx
s
dt (5.8)
The increment in momentum of molecules in t he bo x due t o f orces is
t heref ore a S tr at ono vich s t ochas tic increment, wit h s t ochas tic portion dQ in-
tended t o appro ximate t he micro-integ r ated effect of t he com plicated but
continuous micro-time-scale wiggling of t he individual molecules near t he
surf ace of t he bo x. The use of a S tr at ono vich increment, r at her t han an It ô
increment is due t o t he f act t hat t he r andom increment is a model f or t he de-
tails of a r apidl y fluctuating continuous function, and t he r ules of s tandar d
calculus should appl y t o increments deriv ed from such a process.
Time/dt
Re ff
0 1
Typical Increment
If w e wish t o eliminate t he s t ochas tic portion of t his model w e can tak e
an expectation wit h respect t o t he r andom increment t o find an “effectiv e
a v er ag e f orce” during t he increment. When t he s t ochas tic portion is no t t oo
lar g e, t his will produce a reasonable deter minis tic appro ximate model, in
69
5.4. The Model
addition w e assume t hat density is relativ el y unif or m across t he bar , which
is ho w w e proceed.
N ext, w e ar gue t hat no t onl y is t he r
e
v alue dependent on center of mass
displacement, but also on tem per ature, and t he v elocity difference betw een
t he centers of mass of adjacent bo x es. The ph ysical ar gument f or t his model
goes as f ollo w s:
Suppose t hat w e ha v e tw o bo x es of molecules next t o each o t her and
t hat t he firs t one is ins tantaneousl y heated t o a higher tem per ature t han t he
o t her . The molecules in t his bo x no w ha v e more ener gy and can on a micro-
increment of time occup y a lar g er v olume because t he y can reach higher int o
t he po tential ener gy field produced b y at oms in neighboring bo x es. Thus
molecules near t he boundar y will immediatel y begin t o expand in space,
pushing har der on t he neighboring bo x. The effectiv e dis tance, which de-
scribes t he f orces betw een bo x es, is t heref ore a function of tem per ature. The
f act t hat t his occurs on a micro time increment is due t o t he f act t hat tem-
per ature is a non-local property of t he entire bo x. Similar l y , center of mass
v elocity is a non-local property , an a v er ag e o v er all t he molecules in t he bo x.
The assum p tion is, t heref ore, t hat molecules on t he boundar y betw een t he
bo x es couple t o se v er al molecules wit hin t he v olume and t hese molecules
all couple indirectl y t o t he center of mass ener gy b y t her mal mixing. This
coupling t o center of mass ener gy allo w s t he molecules t o occup y higher
po tential ener gy s tates in t he boundar y region t han t he y o t her wise w ould,
which means t he y can ex ert lar g er f orces, altering t he effectiv e dis tance be-
tw een bo x es on a time scale much f as ter t han t he mo v ement of t he center of
mass.
If w e linear l y T a y lor expand t he r
e
v alues in t he v ariables T; x; v w e can as-
sume t he f ollo wing dimensionless relationship (f or inter action betw een bo x
1 & 2 f or exam ple):
δ r
e
12
x
m
=
A
x
δ= 2
( u
2
u
1
)+
B
v
v
m
( v
2
v
1
)+
C
T
T
0
( T
2
T
1
) (5.9)
In t his expression, δ r
e
12
is t he chang e in r
e
from its equilibrium v alue, x
m
is t he molecular unit lengt h, δ= 2 is t he widt h of t he R VE, and is also t he equi-
librium separ ation betw een successiv e u measurements, v
m
is t he molecular
scale v elocity unit, and T
0
is t he initial tem per ature of t he sys tem. The quan-
tities u
2
and u
1
are t he displacements of t he centers of mass of t he molecules
t hat s tarted in slice 1,2. And A
x
; B
v
; C
T
are dimensionless cons tants. In later
anal ysis w e will define t he model in ter ms of t he bar time scale t
w
which
70
5.4. The Model
will fix t he numerical v alue of t he coefficient A
x
exactl y . Our ques tion t hen
becomes what are B
v
, and C
T
. Thinking firs t about t he tem per ature depen-
dence, if w e con v ert all t he initial center of mass v elocity int o t her mal en-
er gy using numerical v alues typical f or t hese sim ulations and ignoring po-
tential ener gy t her mal capacity , t hereb y o v er -es timating t he dT increment
conser v ativ el y , w e ha v e 3= 2n( T + dT)= T = ( 3= 2n T + n v
2
0
= 2)= T ! dT= T =
2= 3( 0: 1)
2
= 2=( 0: 1) 0: 033 . W e t heref ore conclude t hat chang es in tem per a-
ture are no t im portant f or t he problem at hand as t he tem per ature will be
cons tant t o wit hin a f e w percent t hroughout t he experiment and w e expect
C
T
= O( 1) so w e drop t he T dependence ter m f or t his model.
F inall y , w e tak e B
v
t o be of t he or der
t
m
t
0
where t
0
is t he time it tak es f or a
t her mal wiggle at t he center of mass t o ha v e a probability
1
e
of propag ating
t o t he edg e of t he bo x b y Gaussian diffusion since it is t he coupling of t her -
mal wiggles on t he edg e of t he bo x t o center of mass mo tion t hat causes t he
v elocity dependence of r
e
. If w e assume t he center of mass is near t he center
of t he bo x, t hen t he time it tak es f or a 1D Gaussian r andom w alk t o ha v e 1= e
probability of exiting a bo x of side δ= 2 can be calculated easil y as.
t
e
=
Q( δ= 2)
2
m
m
T t
m
m
m
is t he molecular mass which is tak en as 1 in our simulation, t
m
is t he
molecular mean free pat h time which w e tak e t o be of t he or der of t he molec-
ular time unit. The cons tant Q =
1
8
[erf
1
( 1 1= e)]
2
0: 308 is due t o t he 1= e
probability chosen f or exiting t he bo x.
Ho w e v er , no tice t hat due t o t he non-zero r ang e of t he LJ f orces when
t he bo x is small enough, t her mal mo tion t hroughout t he bo x immediatel y
affects t he surf ace. W e t heref ore propose t he f ollo wing modification which
has cons tant beha vior asym p t o ticall y f or small bo x es and similar asym p t o tic
beha vior f or lar g e bo x es as t he pre vious result:
t
0
=
Q( 1+( δ= 2= x
0
)
2
) x
2
0
m
m
2t
m
T
Where x
0
is of t he or der of t he bo x size at which t he Lennar d-Jones f orces
per meate t he entire bo x.
Subs tituting t he expression f or t he chang e in r
e
(equation 5.9 ) int o t he
model f or dP (equation 5.8 ) and dividing b y dt w e ha v e:
71
5.5. Equiv alent Continuum Model
dP
dt
= ( ρ A dx
s
) K x
m
(
A
x
( u
2
2u
1
+ u
0
)
δ= 2
+
B
v
v
m
2t
2
m
T
m
m
Q( 1+(
δ
2x
0
)
2
) x
2
0
( v
2
2v
1
+ v
0
)
)
(5.10)
Here K 57: 15 is t he s tiffness f or t he quadr atic T a y lor expansion of t he
Lennar d-Jones po tential about t he equilibrium location at its minimum, in
t he dimensionless LJ “units. ” W e let dx
s
= x
0
in our s tatis tical model . If
in defining t he model abo v e, w e express time in units of t
w
and dis tance
in units of δ= 2 t hen necessaril y t he square of t he w a v e v elocity ρ A x
0
K A
x
=
100 . W e specified t his f or m in our s tatis tical model, and proceeded t o fit
t he beha vior of t he model t o t he obser v ed w a v e deca y time-series, including
inf erring bo t h B
v
and t he precise v alue of t
w
f or each bar .
5.5 Eq uiv alent Continuum Model
T o mak e t he tr ansition from a discrete model f or a small fix ed set of obser v a-
tions in a tin y bar (where t he R VE w as an appreciable fr action of t he lengt h
on t he or der of 1/20) t o a continuum model f or an unkno wn number of mea-
surements o v er a v er y long dis tance in a v er y t hin bar , w e can reinter pret t he
alg ebr aic quantities in equation 5.10 b y allo wing δ
= 0 when measured at t he
bulk scale. W e begin b y re writing t he dimensionless model in equation 5.10
using ne w dimensionless v ariables in ter ms of fr actions of t he bulk scale dis-
tance, time, and mass. Also w e f act or t he mass of t he R VE slice out of t he
momentum v ariable P = dm
′
v .
dm
′
dv
dt
= ( ρ A m
m
dx
s
) K x
m
(
A
x
( u
2
2u
1
+ u
0
)
δ= 2
+
B
v
2t
2
m
T
m
m
Q( 1+(
δ
2x
0
)
2
) x
2
0
( v
2
2v
1
+ v
0
)
v
m
)
(5.11)
W e begin b y no ticing t hat dm
′
= ρ A m
m
δ= 2 t he dimensionless mass of t he
slice, and dividing bo t h sides b y t his quantity t o g et a model t hat is in v ariant
t o t he mass scale:
72
5.5. Equiv alent Continuum Model
dv
dt
= dx
s
K x
m
(
A
x
( u
2
2u
1
+ u
0
)
( δ= 2)
2
+
B
v
2t
2
m
T
m
m
Q( 1+(
δ
2x
0
)
2
) x
2
0
( v
2
2v
1
+ v
0
)
v
m
( δ= 2)
2
δ= 2
)
(5.12)
N ext w e con v ert t he entire equation t o dimensional f or m on t he bulk
scale, b y multipl ying all lengt hs, and times b y dimensional quantities dL =
ε
L
L and dt = ε
t
t
. The coefficient K is multiplied b y 1= dt
2
as it expresses
a f orce per unit mass per unit displacement as a fr action of t he molecular
units. On t he bulk scale x
m
= 1dL
dL
dt
2
dv
dt
= dL dx
s
dL
K
dt
2
(
A
x
( u
2
2u
1
+ u
0
)
( δ= 2)
2
dL
dL
2
+
B
v
t
m
t
0
( T; δ)
( v
2
2v
1
+ v
0
)
v
m
( δ= 2)
2
1
dL
2
δ= 2dL
)
(5.13)
N o w w e inter pret dL dx = dL δ= 2 as an infinitesimal displacement at t he
bulk scale, and re write t he r atios of t he differences using t he definition of
t he deriv ativ e in IS T dy= dx
= ( y
1
y
0
)= dx . W e also cancel t he f act or dL= dt
2
.
dv
dt
= dx
s
dL K
(
A
x
@
2
u
@ x
2
1
dL
+ B
v
t
m
t
0
( T; δ)
@
2
v
@ x
2
δ
2dL
)
(5.14)
canceling t he dL w e arriv e at t he continuum model:
dv
dt
= dx
s
K
(
A
x
@
2
u
@ x
2
+
B
v
δ
2
t
m
t
0
( T; δ)
@
2
v
@ x
2
)
(5.15)
F inall y , assuming t he lengt h scale is fix ed at L , t he bar scale lengt h, w e can
rescale ag ain t he time scale, replacing t he dimensionless t wit h ε
t 2
t
′
and t he
dimensionless v b y v
′
= ε
t 2
, and multipl ying bo t h sides b y ε
2
t 2
wit h ε
t 2
chosen
so t hat dx
s
K A
x
ε
2
t 2
= 1 . The model in sim plified bar lengt h and time scales is:
dv
dt
′
=
@
2
u
@ x
2
+
B
v
ε
t
dx
s
K δ
2
t
m
t
0
( T; δ)
@
2
v
′
@ x
2
(5.16)
N o w w e examine t he coefficient B
v
ε
t 2
dx
s
K δ= 2 , and ask whet her t his is
infinitesimal, appreciable, or unlimited? Clear l y K and dx
s
are limited as
73
5.5. Equiv alent Continuum Model
K 57: 15 and dx
s
= ord( 1) w as t he multiplier f or t he molecular lengt h
t hat g a v e t he surf ace region t hickness of t he R VE. S tatis ticall y B
v
w as O( 1)
(see chap ter 6 ) and δ= 2 > 1 is an appreciable finite R VE t hickness as a fr ac-
tion of t he molecular widt h. The ques tion remains onl y ho w lar g e is ε
t 2
=
1=( dx
s
K A
x
) but since w e expect A
x
= ord( 1) it appears t hat t his whole coef-
ficient should be considered near -s tandar d. Incor por ating t he cons tants int o
a single cons tant t hat depends on t he o bser v ational scale δ , B
′
v
( δ) w e ha v e
t he final model:
@ v
@ t
=
@
2
u
@ x
2
+ B
′
v
( δ)
t
m
t
0
( T; δ)
@
2
v
@ x
2
(5.17)
It should be no ted t hat t his model depends on δ which is a measure of
t he r atio of t he scale of obser v ations t o t he scale of t he inter -molecular spac-
ing. When δ is lar g e t he asym p t o tic beha vior of t he coefficient in front of
@
2
v
@ x
2
is proportional t o 1= δ , so t hat at lar g e obser v ational sizes t he momentum
diffusion effect becomes negligible. F or seismic w a v es tr a v eling in t he cr us t
of t he eart h, t her mall y driv en momentum diffusion is com pletel y ignor able.
Ho w e v er f or w a v es tr a v eling t hrough a small cr ys tal and arriving at a tin y
micro-electro-mechanical (MEMS) accelerometer t he effect could be im por -
tant. This kno w ledg e ma y help inf or m sensor based cr ack f or mation and
cr ack propag ation models in s tr uctur al monit oring or similar applications.
This continuum equation came about b y reinter preting an inter mediate
scale molecular model, which w e de v eloped in t he context of a v er y w ell de-
fined MD simulation, where t he obser v ational scale w as kno wn exactl y . This
continuum equation can be used when w e ha v e a t hin molecular bar which
is much long er t han its t hickness, and where w e ha v e a measurement ap-
par atus which measures displacements or v elocities o v er a lengt h of t he bar
which is small relativ e t o t he full bar lengt h. In such a context, as described
in chap ter 3 it is reasonable t o treat t he measurement region as an infinitesi-
mal region relativ e t o t he t o tal lengt h, and t o appl y t he IS T tech niques when
building our model. This allo w s us t o tr ansition from discrete difference
equations f or a small set of measurements t o a differential equation, but onl y
so long as t he inter mediate asym p t o tics t hat are assumed b y t hat tr ansition
hold t o sufficient appro ximation. The IS T based approach mak es t his type of
reasoning s tr aightf or w ar d, replacing confusing limiting ideas about a non-
ph ysical continuum idealization, int o a combination of alg ebr aic mat hemat-
ics and asym p t o tic appro ximate reasoning about real ph ysical quantities.
74
5.6. Conclusion
5.6 Conclusion
Our goal w as t o g ener ate a mat hematical model t hat explains t he dissipation
of w a v e ener gy in a sys tem t hat is decoupled from t he en vironment, where
inter nal entrop y g ener ation is t he onl y mechanism of dissipation. By using
t he mat hematical t ool of N ons tandar d Anal ysis as an aid t o our f or mulation,
w e w ere able t o account f or and describe t he separ ation of scales t hat occurs
in a sys tem where molecules under go detailed undulations t hat are much
f as ter t han t he chang es in t he regionall y a v er ag ed quantities deriv ed from
t hose molecules. T o couple t he detailed dynamics t o t he a v er ag ed dynamics
w e made ref erence t o ener g etic ar guments t hat pro vided a scaling la w v alid
asym p t o ticall y f or lar g e regions. By modification of t his asym p t o tic la w t o
ha v e cons tant beha vior f or smaller regions w e arriv ed at an equation which,
when fit t o data, pro vides a quite accur ate descrip tion of t he r ate of dissi-
pation of w a v es across obser v ational scales from betw een about 1000 and
100000 molecules, and a wide r ang e of tem per atures in t he solid r ang e f or
Lennar d Jonesium. In t he next chap ter , I will discuss t he fit of t his model t o
t he actual data from MD simulations, and t he Ba y esian s tatis tical techniques
used t o com pare t he model t o t he data.
75
Chapter 6
Comparing Dissipating W av e M odel to
D etailed Timeseries
T o under s t and God’ s t houghts w e mus t s tudy s t atis tics, f or t hese ar e t he
measur e of his pur pose.
— Florence Nighting ale (attr .)
In chap ter 5 I de v eloped a model f or bulk -scale kinetic ener gy dissipation in
a vibr ating molecular bar . The model has t hree principal deg rees of freedom
t hat should be considered in or der t o fit t he data. The mos t im portant and
fundamental deg ree of freedom is B
v
t he multiplier f or t he momentum dis-
sipation process. The assum p tions of t he model im pl y t hat t his should be a
number in t he vicinity of 1 t hank s t o t he dimensional anal ysis t hat w ent int o
t he model cons tr uction. The o t her deg rees of freedom f or fit-t o-data are t he
precise time-scale t
w
which depends on t he precise size of t he bar and t he
w a v e speed, and t he t her mal kinetic noise le v el K T which is a small multiple
of 20T= 2 since w e are measuring 20 deg rees of freedom each of which con-
tributes T= 2 t o t he kinetic ener gy on a v er ag e. In addition t o t hese per -bar
deg rees of freedom, t he o v er all model has a lengt h scale x
0
which specifies
t he asym p t o tic baha vior of t he model f or small bars whose obser v ational
slices are O( x
0
) wide or smaller .
Using Mar k o v Chain Monte Car lo met hods, w e fit a Ba y esian s tatis tical
model f or t he error betw een t he predictions from t he ODE model abo v e and
t he actual s tatis tics of t he molecular dynamics time-series, giving v alues f or
t he B
v
, t
w
, K T , and x
0
par ameters. This model w as fit t o 40 r andoml y chosen
time-series from t he t o tal lis t of a v ailable experiments. The model had x
0
as
76
a global par ameter , and B
v
, t
w
, K T , and E
0
as bar specific par ameters. The pa-
r ameter K T is t he minimum t her mall y deriv ed ener gy t hat w ould be present
e v en if t he w a v e w ere absent and w as specified as a multiple of 20T= 2 . The
par ameter E
0
is a fitting par ameter t hat multiplies t he initial w a v e ener gy t o
account f or noise in t he initial conditions of t he molecular dynamics simula-
tion t hat result in small differences betw een t he actual and nominal v alues.
A prior w as assigned f or each of t hese par ameters which w as common t o t he
v alues across all bars. The B
v
v alues w ere giv en a common log-nor mal prior ,
and t he s tandar d de viation of t he under l ying n or mal w as a h yper par ameter
and giv en its o wn exponential prior wit h mean v aluelog( 1: 5) .
In figure 6 (t op) w e see t hat t he mean B
v
v alues f or t he 40 simulations
f or m an ensemble whose dis tribution is significantl y narro w er t han t he prior
used, such t hat B
v
2 [ 0: 5; 2] f or t he ma jority of cases. V ariability in B
v
ma y
in part be accounted f or b y t he details of t he size and arr ang ement of pol y -
cr ys talline domains in t he bar . In t he lo w er g r aph w e see t he dis tribution of
pos terior sam ples of B
v
f or f our exam ple bars, b y color , indicating t hat no t
onl y are t he mean v alues f or each bar concentr ated in t he r ang e [ 0: 5; 2] but
t he individual bars also ha v e individual sam ples t hat are concentr ated in a
subset of t his r ang e.
Using t he pos terior means of t he par ameters, w e g r aph t he w a v e deca y
time-series and t he v elocity time-series f or t he t hir d bar slice in f our r an-
doml y chosen bars. Clear l y t he beha vior is qualitativ el y and quantitativ el y
correct. Alt hough t he model is fit t o t he t o tal w a v e deca y ener gy time-series
onl y , t he equations also correctl y predict t he details of t he f orces on slice
3, resulting in t he proper v elocity time-series sho wn on t he right. The cor -
rectness of t he v elocity time-series giv es us confidence t hat t he momentum
diffusion process w e model is subs tantiall y correct.
As can be seen in figure 6.2 in our f our r andoml y chosen exam ple bars,
when t he pos terior mean v alues are used f or B
v
and t
w
no t onl y is t he pre-
dicted time-series of kinetic ener gy closel y matched t o t he obser v ed data, but
t he time-series f or t he predicted v elocity in t he t hir d spatial measurement is
closel y matched t o t he obser v ed v alue.
F inall y , w e examine trends in t he mean B
v
v alues wit h eit her T or δ and
find t hat no s trong trends exis t.
77
Prior
Posterior
0
1
2
3
0.25 0.50 1.00 2.00 4.00
Bv (log scale)
density
Distribution of mean Bv values
across 40 bars
0
1
2
3
0.25 0.50 1.00 2.00 4.00
Bv (log scale)
density/4
Distribution of Bv samples
across 4 example bars
F igure 6.1: Dis tribution of t he pos terior mean of t he B
v
v alues f or each of
40 exam ple bars, com pared t o t he prior dis tribution (t op). And dis tribution
of t he B
v
sam ples f or t he same 4 exam ple bars as sho wn belo w (individual
colors, bo tt om). Each bar has its o wn B
v
but t he y are concentr ated in t he re-
gion [0.5,2]. The B
v
v alue dis tribution in an y giv en bar is probabl y influenced
b y t he r andom pol y cr ys talline s tr ucture induced b y t he initialization condi-
tions. The prior is e v aluated at t he pos terior mean of t he h yper par ameters
t hat define it.
78
Predicted
Measured
−6
−4
−2
0
0 5 10 15
time/t_w
log(Wave Energy / Initial)
Wave Energy
Predicted vs Measured
−0.2
−0.1
0.0
0.1
0.2
0 5 10 15
time/t_w
V[3] measured or predicted
Comparision of V[3]
Predicted vs Measured
Predicted
Measured
−6
−4
−2
0
0 5 10 15
time/t_w
log(Wave Energy / Initial)
Wave Energy
Predicted vs Measured
−0.2
−0.1
0.0
0.1
0.2
0 5 10 15
time/t_w
V[3] measured or predicted
Comparision of V[3]
Predicted vs Measured
Predicted
Measured
−6
−4
−2
0
0 5 10 15
time/t_w
log(Wave Energy / Initial)
Wave Energy
Predicted vs Measured
−0.2
−0.1
0.0
0.1
0.2
0 5 10 15
time/t_w
V[3] measured or predicted
Comparision of V[3]
Predicted vs Measured
Predicted
Measured
−6
−4
−2
0
0 5 10 15
time/t_w
log(Wave Energy / Initial)
Wave Energy
Predicted vs Measured
−0.2
−0.1
0.0
0.1
0.2
0 5 10 15
time/t_w
V[3] measured or predicted
Comparision of V[3]
Predicted vs Measured
F igure 6.2: Molecular Dynamics v s t he 20 D model in 4 r andoml y chosen
exam ples.
79
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.0
0.5
1.0
1.5
2.0
0.0 0.1 0.2
T
Bv
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.0
0.5
1.0
1.5
2.0
0 10 20 30
x
Bv
F igure 6.3: Examination of residual trends in B
v
re v eal no sys tematic biases
wit h v ar ying T or δ across t he 40 bars used, indicating t hat t he a v er ag e effects
of tem per ature and size are w ell accounted f or in t he model.
80
6.1. Issues W it h F itting Ba y esian Models t o Dynamics
6.1 Issues Wit h Fitting Ba y esian Models to
Dynamics
The Chosen T echniq ue
In or der t o deter mine t he B
v
v alues in t he model f or w a v e dissipation, it w as
necessar y t o r un t he ODE predictions and tr y t o disco v er which B
v
v alues
w ere consis tent wit h actual beha vior , and t hen whet her t hose v alues w ere
consis tent wit h t he assum p tions of t he dimensional anal ysis and modeling
t hat led t o t he model.
In attem p ting t o fit t he Ba y esian s tatis tical model certain interes ting as-
pects of t he problem became clear and w ere considered illuminating f or t he
nature of fitting dynamic models in g ener al. The Ba y esian model w as fit us-
ing Mar k o v Chain Monte Car lo (MCMC) and had t he f ollo wing s tr ucture.
The ODE w as r un t o produce position and v elocity es timates at all of t he
obser v ed time points f or all 40 exam ple bars. The predicted w a v e ener gy
w as calculated as t he sum W
p
=
∑
20
i= 1
v
i
( t)
2
= 2 and t his ener gy divided b y
its initial v alue t o g et an ener gy on t he scale of t he initial conditions. The
mass v alues f or each deg ree of freedom w ere assumed t o be t he same, and
are nor malized a w a y b y t his rescaling.
The w a v e ener gy data from MD simulations w as similar l y rescaled as a
fr action of t he nominal initial e ner gy W
d
= ( W
i
=( 6= 20N v
2
init
= 2))
i
f or all of
t he bars index ed b y i .
A prior f or bar timescale t
i
w as specified f or each bar using a nominal
v alue f or t
w
= 5( N= 10)
( 1= 3)
= v
w
as t he center of t he prior dis tribution wit h
coefficient of v ariation of 0: 133 , v
w
w as tak en t o be t he nominal v alue
√
K
ρ
7: 38 .
The lik elihood f or a giv en set of t
i
and B
v
i
w as specified as independent
Gaussian processes on tr ansf or med errors in t he log r atio of t he w a v e ener -
gies.
p( q
N
( p
t
(log(( W
p
+ K T)=( E
i
W
d
))= S; n))j B
v
; t
i
; S;:::) = N( 0; C) (6.1)
where p
t
is t he CDF f or t he s tandar d s tudent t dis tribution ( wit h n de-
g rees of freedom, which w as giv en a h yper -prior), q
N
is t he quantile func-
tion of t he s tandar d nor mal (in v erse of t he CDF), S is an error scale which
includes effects bo t h f or modeling error , and t her mal noise, E
i
is an initial en-
er gy multiplier accounting f or differences in t he initial w a v e ener gy from its
81
6.1. Issues W it h F itting Ba y esian Models t o Dynamics
nominal v alue, and K T is a t her mal base ener gy le v el added t o t he noise-free
predict or .
This might be called a “Gaussian Copula t-Process” since t he mar ginal
dis tribution at an y time point of t he log r atio is t dis tributed, but t he depen-
dency in t he model comes from a Gaussian process t hrough t he probability
integ r al tr ansf or m. The s tudent t tr ansf or mation w as used in or der t o ac-
count f or t he f act t hat as t he kinetic ener gy reaches a minimum at each cy -
cle (when t he w a v es o v er lap in t he center of t he bar) t he log r atio becomes
log( W
p
+ K T)log( E
i
W
d
) and W
d
is a small positiv e noisy v alue, and t here-
f ore outliers are t o be expected and should no t ex cessiv el y affect t he fit.
The co v ariance matrix C w as built using t he f ollo wing co v ariance func-
tion:
cov( t
1
; t
2
) = 0: 1exp((( t
1
t
2
)= 0: 2)
2
)+ 0: 9δ( t
1
; t
2
)sin
2
( π t
1
= 2) (6.2)
This co v ariance function reflects our kno w ledg e t hat a w ell-fitting model
will predict t he peak s of t he w a v e ener gy v er y w ell at time points 0; 2; 4:::
when t hesin ter m is at a minimum so t hat t he mar ginal s tandar d de viation
is also at a minimum at t hose time points, and will ha v e lar g er v ariation a w a y
from t hese time points, wit h maximum v ariation near t he ener gy minima.
F urt her more it predicts correlations betw een errors t hat are nearb y in time
on t he or der of 0.2 times t he cy cle lengt h. The o v er all scale of t he v ariation
is already accounted f or b y S in equation 6.1 .
F itting t he MCMC model required tuning t his co v ariance function, and
t he jum ping scales f or t he metropolis diffusion. The “mcmc” packag e (Ge y er
and Johnson 2013 ) w as used in t he s tatis tical sof tw are R t o im plement t he
Mar k o v Chain because it allo w ed f or defining t he log pos terior density via
r unning t he ODE sol v er t o calculate t he errors. Because of t he e v aluation of
40 ODE based time-series twice f or each jum p (f or w ar d and backw ar ds t o
ensure detailed balance) t he com putation proceeded slo w l y on t he deskt op
com puter on which it w as r un, ultimatel y taking about a w eek f or t he final
r un t o g ener ate on t he or der of 10,000 sam ples.
In carr ying out t his s tatis tical in v es tig ation, se v er al t hings became clear
t hrough experience. The firs t w as t hat in man y problems, especiall y prob-
lems in dynamics, t here is no a-priori correct lik elihood. In t hese cases, t he
procedure f or fitting a s tatis tical model ha v e a different meaning t han t he y
do in a more w ell defined context such as a r andomized controlled clinical
82
6.1. Issues W it h F itting Ba y esian Models t o Dynamics
trial of one condition v s ano t her . W e choose a model t o fit some pur pose,
r at her t han because of some particular prior belief about t he f or m of t he
model. F urt her more, when t he model is an indirect or appro ximate model of
some more detailed phenomenon, as wit h t he dissipativ e w a v e model of t he
more detailed MD simulation, t here do no t exis t “tr ue” v alues of t he par am-
eters, especiall y B
v
which is t he main non-nuisance par ameter . As pointed
out in Kilmins ter and Machete ( 2005 ) under t his condition, it mak es sense
t o adjus t t he s tatis tical process t o achie v e some goal. Alt hough t he y specify
some particular goals rele v ant t o t heir problem, in g ener al w e will ha v e some
types of fits considered accep table f or our pur poses, and o t her types of fits
considered unaccep table.
The criteria f or choosing t he particular co v ariance function here w ere t hat
t he s tatis tical model should no t allo w or pref er f or t he timescales t o become
significantl y mismatched betw een t he data and t he predictiv e model, and
t hat t he peak w a v e ener gy v alues should be fit w ell f or all of t he a v ailable
obser v ation time, so t hat t he deca y r ate of t he w a v e ener gy w as w ell pre-
dicted. The f act t hat cons tr aining onl y t he w a v e ener gy deca y led also t o
subs tantiall y correct time-series f or t he individual v elocities at each point
on t he bar lends credence t o t he f or m of t he model and t he exi s tence of t he
ph ysical process t hat w as considered.
R elationship to Appro ximate Ba y esian Com put ation (ABC)
It should be mentioned t hat t here is a sense in which t his procedure fits
wit hin t he Appro ximate Ba y esian Com putation (ABC) fr ame w or k (f or a re-
vie w see eg. Marin et al. ( 2012 )). Each bar has 20 v elocity measurements at
each of about 80 time p oints bo t h in t he data and in t he ODE model. There
are also displacement measurement in t he ODE model, but t hese are no t ob-
ser v ed in t he MD model wit h fix ed Eulerian measurement bo x es. W e ha v e
no w a y t o write do wn t he lik elihood of particular data giv en B
v
and o t her
nuisance par ameters since t here is no a-priori r andom data g ener ating pro-
cess t o pref er , but w e can simulate pseudo-data from t he appro ximate model
so t hat w e g et t he full 20 dimensional v elocity v s. time tr a ject ories under an
assumed v alue of B
v
and suitable initial conditions. As is common in ABC,
r at her t han fitting our simulation t o t he entire data set of obser v ation, in-
cluding t he entire 20 D v elocity time-series, w e choose a summar y s tatis tic
equal t o t he time-series of t o tal w a v e kinetic ener gy . Ins tead of a single ε
v alue giving a t oler ance f or de viation betw een a small dimensional sum-
83
6.1. Issues W it h F itting Ba y esian Models t o Dynamics
mar y s tatis tic of actual data (from MD) and predicted pseudo-data (from
t he ODE model) w e specify a kind of k er nel in t his 80 dimensional w a v e
ener gy time-series using a particular tr ansf or mation of a Gaussian process
wit h a particular co v ariance s tr ucture. In doing so, w e cons tr uct a kind of
k er nel-density -lik elihood in which w e ha v e im plicitl y described t he char ac-
ter of t he modeling errors t hat our ODE model is allo w ed t o introduce. The
result leads t o ex cellent inf erence on B
v
and good ag reement wit h t he w a v e
ener gy deca y time-series as w ell as t he time-series in t he 20 additional v eloc-
ity dimensions f or each bar .
Efficiency of MCMC Calculation wit h Expensiv e Simulations
Ano t her is sue which became ob vious in r unning t he MCMC simulation w as
t hat t he log pos terior density (LPD) w as calculated (up t o a nor malizing con-
s tant offset) at a lar g e number of points, but t hese v alues w ere t hro wn a w a y
immediatel y af ter each jum p. Each e v aluation in v ol v ed calculating t he full
time ser ies of 40 different ODE models and t hen com paring t hem t o t he
data. The results are continuousl y dependent on e v er y par ameter , so t hat
tw o e v aluations which are close in par ameter space mus t be close in LPD.
By s t oring t he LPD v alue wit h each e v aluation, w e can quickl y appro ximate
t he LPD when a proposal brings us close t o a region w e ha v e pre viousl y seen
(alt hough t his w as no t done in our fit). V arious proposals f or met hods re-
lated t o t his insight ha v e been made in t he liter ature. F or exam ple in T addy ,
Lee, and Sansó ( 2009 ) t he aut hors sugg es t a met hod whereb y t he support of
t he par ameter space is discretized ont o a set of v alues where t he simulation
results are kno wn, and t hen t he pos terior dis tribution o v er t hose discrete
locations is deter mined using resam pling. In Higdon et al. ( 2004 ) t he au-
t hors propose fitting a surrog ate s tatis tical model, so t hat pseudo-data which
w ould nor mall y come from expensiv e simulations ins tead comes from some
f or m of inter polation of a set of fix ed simulation r uns. In essence t heir inter -
polation technique in v ol v es t he use of a Gaussian process prior on t he un-
obser v ed simulation output conditional on t he obser v ed simulation output,
t his can be vie w ed as a f or m of inter polation. A useful discussion of Gaus-
sian process models, and t heir equiv alence t o v arious f or ms of par ametric
spline or r adial basis function smoo t hing techniques is a v ailable in MacKa y
( 2003 ). In certain circums tances it can be useful t o model unkno wn func-
tions as nonlinear tr ansf or mations of Gaussian processes, so called “w ar ped
Gaussian processes” (Snelson, Rasmussen, and Ghahr amani 2004 ). This can
84
6.1. Issues W it h F itting Ba y esian Models t o Dynamics
help t o account f or non-Gaussian error dis tributions, and is a g ener alization
of t he technique used t o account f or outliers when building t he error model
in t he dissipativ e bar experiment.
85
Chapter 7
F uture Directions f or Continuum
M odeling of Sand Grains During
Liquef action
W e ha v e no right t o assume t hat an y phy sical la w s e xis t, or if t he y ha v e
e xis t ed up until now , t hat t he y will continue t o e xis t in a similar manner
in t he futur e.
— Max Planck The U niv er se in t he Light of Modern Phy sics (1931)
The model de v eloped in chap ter 5 w as intended t o be a tes tbed f or modeling
ideas de v eloped in chap ter 3 and f or fitting of Ba y esian s tatis tical models t o
dynamics. The original mo tiv ation f or de v eloping t he modeling techniques
ho w e v er , w as t o explore ho w sand g r ains could be modeled in a continuum
fr ame w or k in a consis tent w a y so as t o couple a g r ain mo tion continuum
model t o t he pore-pressure equations of chap ter 4 . Clear l y , g r ains of sand
are discrete, y et w e ha v e v er y limited inf or mation about t heir individual size,
shape, and arr ang ement so t hat models in v ol ving individual discrete parti-
cles are artificiall y sim plified t o spheres, incapable of dealing wit h realis tic
shapes, or prohibitiv el y com putationall y expensiv e. Also, w e ha v e usuall y
at bes t coarse measurements of displacements corresponding t o spatial a v -
er ag es o v er man y man y g r ains so t hat t he details of g r ain mo tion are g ener -
all y non-f alsifiable an yw a y . In t his setting, when inter preting a continuum
model as spatiall y a v er ag ed o v er spatial scales much smaller t han t he full
scale, a continuum model seems ideal so long as w e can de v elop one which
86
7.1. The Setting
giv es accur ate predictions while maintaining a causal connection t o t he sub-
scale processes t hat produce t he a v er ag es.
In t his chap ter , w e consider some mechanical issues t hat should be tak en
int o account when tr ying t o build such a continuum model f or g r ain mo tion.
7.1 The Setting
Imagine a slice of soil similar t o t he one used t o de v elop t he w ater flo w model
in chap ter 4 . The cross sectional area is much lar g er t han t he cross section of
a single g r ain, and t he t hickness is much smaller t han t he o v er all t hickness
of t he soil deposit, but t he t hickness is modeled as some s tandar d or non-
s tandar d multiple of t he typical sand g r ain diameter D , w e will in v es tig ate
bo t h op tions. The soil slice is in dynamic equilibrium so t hat t he net f orce on
t he slice is zero. The mass of t he slice is M = A ρ
g
( 1 φ) dz+ A ρ
w
sφ dz+ ε
m
where s is t he deg ree of satur ation wit h w ater (air is neglected) and ε
m
ac-
counts f or t he errors induced b y small fluctuations, f or exam ple when g r ains
cross t he slice boundar y . Because of t he lar g e cross sectional extent t here
will be a lar g e absolute number of g r ains o v er which w e a v er ag e so t hat t he
centr al limit t heorem im plies an appro ximatel y nor mall y dis tributed error
in t he a v er ag e ε
m
N( 0; σ) wit h σ small relativ e t o t he o v er all M v alue. If
t he slice a v er ag es o v er a nons tandar d number of g r ains t hen w e expect σ= M
is infinitesimal, whereas f or a s tandar d multiple t he scale of t he s tatis tical
noise will be appreciable.
W e d ivide t he slice int o a possibl y nons tandar d number of sub-slices.
There are a certain number of sub-slices near t he bo tt om and t he t op which
define a v olume “at t he surf ace” of t he slice. The inclination at t his point is t o
treat t he “surf ace v olume” as infinitesimal relativ e t o t he t o tal slice v olume,
but t his assum p tion ma y tur n out t o be untenable. The reason f or defin-
ing a “near -surf ace v olume” r at her t han sim pl y a “surf ace area” is t hat all
ph ysical objects occup y some v olume. The contacts betw een slices are no t
all planar along a surf ace area, and it is t he contact f orces which will driv e
t he mechanical beha vior . See figure 7.1 f or a g r aphical representation.
The net momentum of t he g r ains in t he slice of interes t is affected b y
f orces tr ansmitted from g r ains in t he near -surf ace v olume of neighboring
slices t o g r ains in t he near -surf ace v olume of t he slice of interes t via contact
f orces. These surf ace f orces are t hen tr ansmitted int o t he v olume where t he y
cause t he acceler ation of t he whole slice. If ins tead w e model t he slices as
87
7.1. The Setting
dz
ε dz
An Infinitesimal slice of soil with
Infinitesimal sub-slices
F igure 7.1: An illus tr ation of t he basic concep t of an infinitesimal slice of
soil of t hickness dz brok en int o infinitesimal sub-slices of t hickness ε dz . The
“near surf ace v olume” w ould be per haps one or tw o sub-slices t hick. (ac-
tual scales relativ e t o g r ain size chosen f or con v enience of illus tr ation no t
quantitativ e accur acy)
appro ximatel y t he t hickness of a f e w g r ains, so t hat all of t he v olume of t he
slice is “near -surf ace” w e ma y ha v e a more difficult time wit h t he model
since some inter actions extend o v er 5 t o 10 g r ain lengt hs in realis tic g r anular
materials possibl y cutting t hrough entire slices (Ma jmudar and Behring er
2005 ).
N ext w e consider a time dt which is infinitesimal wit h respect t o t he scale
of t he o v er all e v ent of i nteres t. T ypical eart hquak es las t tens t o up t o hun-
dreds of seconds so per haps 1 unit of non-dimensional time representing 30s
is t he or der of magnitude of a reasonable time scale. Deposits of interes t f or
soil liquef action are on t he or der of 10 m deep, and im portant deposits can be
as t hin as 1 cm or less. T reating dx = 0: 01m as infinitesimal relativ e t o 10m
deposit dep t h, sugg es ts t hat t he infinitesimal time should be on t he or der of
dx= v where v is t he bulk scale w a v e speed in sand. T ypical w a v e speeds in
solid miner als of t he same type as sand g r ains are on t he or der of 3000m=s ,
but t he em pirical bulk w a v e speed t hrough loose sand is much lo w er t han
88
7.1. The Setting
3000m=s , and is ma ybe typified b y on t he or der of 300m=s ins tead. The “in-
finitesimal time” is t heref ore on t he or der of : 01= 300 = 3: 3 10
5
s , or 33
microseconds.
Ho w e v er , in 33 microseconds, wit hin t he body of a single sand g r ain, de-
f or mations should tr a v el at t he bulk w a v e speed in t he miner al. The im plied
dis tance is 3000 3: 3 10
5
: 1m or about 100 times t he diameter of a 1mm
sand g r ain. The im plication is t hat w a v es entering a single sand g r ain mus t
be reflected, scattered, and tr ansmitted in such a w a y t hat t he o v er all bulk
w a v e is in essence t he interf erence patter n caused b y all of t his reflection and
scattering.
It’ s kno wn t hat in sand g r ains, f orce chains f or m where se v er al sand
g r ains in a ro w are coupled t og et her tightl y b y s trong contact f orces, and
o t her g r ains are near l y un-loaded (Ma jmudar and Behring er 2005 ). The flux
of momentum t hrough t he boundar y of t he soil slice consis ts of a set of
im pulses produced b y t h e contact f orces in our infinitesimal surf ace v ol-
ume. Along some of t hese s trongl y coupled chains, t he im pulse should
tr a v el quickl y , closer t o t he miner al w a v e speed t han t o t he bulk speed. And,
along o t her less tightl y coupled chains t he im pulses will tr a v el more slo w l y .
Similar l y f or t he net angular momentum tr ansport caused b y t he t or ques due
t o t he off-center nature of t he contact f orces. If a f orce chain is tightl y sup-
ported, it ma y tr ansmit t or que, whereas a t or que on a less tightl y coupled
chain ma y cause chain collapse.
W e consider a coor dinate sys tem where z points upw ar ds t o w ar ds t he
surf ace, x is t o t he right, and y goes a w a y from t he vie w er . W e assume a
v erticall y tr a v eling plane w a v e whose direction of vibr ation is tr ansv erse t o
y in t he x direction.
M dv
x
=
N
b
∑
i= 1
F
x i
dt
N
t
∑
i= 1
F
x i
dt (7.1)
M dv
z
=
N
b
∑
i= 1
F
zi
dt
N
t
∑
i= 1
F
zi
dt M g dt (7.2)
Where in t he abo v e equations, N
b
is t he number of contact f orces at t he
bo tt om surf ace v olume of t he slice and N
t
is t he number of f orces at t he t op
surf ace v olume.
W e assume t hat t he net t or que on t he slice is zero, since t he soil does no t
ro tate as a whole, ho w e v er individual t or ques are definitel y non-zero and
89
7.2. A sliding t hought experiment
can cause individual g r ains t o ro tate which in tur n causes t he ro tation of
o t her g r ains, and ultimatel y t he rearr ang ement of t he particles int o denser
or looser packing.
7.2 A sliding t hought exper iment
Imagine t hen, tw o slices of soil s tack ed on t op of each o t her , wit h t he bo t-
t om slice of soil acceler ating at an im posed horizontal acceler ation a
l
( t) = t
linear in t he time. Let a be non-dimensional in fr actions of t he s teady s tate
acceler ation of t he upper slice, and t non-dimensional in fr actions of t
y
t he
time it tak es f or t he bo tt om slice t o reach t he same acceler ation as t he upper
slice’ s e v entual s teady s tate. In or der f or t his par ameterization t o w or k t here
mus t be a s teady s tate acceler ation f or t he upper slice. Wh y w ould t his be
so?
Ins tead of tw o slices of soil, imagine a con v e y or belt, and a block of w ood
wit h sand paper wr apped around it. The con v e y er belt begins t o acceler ate,
and initiall y it tak es t he block of w ood wit h it. Ev entuall y , as t he con v e y er
belt is acceler ated f as ter and f as ter , t he f orce necessar y t o k eep t he block of
w ood acceler ating wit h t he con v e y or belt ex ceeds what w e t hink of as t he
f orce due t o t he s tatic coefficient of friction, and t he f orce on t he block be-
comes near l y cons tant so t hat t he block acceler a tes at a cons tant asym p t o tic
acceler ation. W e assume a similar t hing will happen wit h our tw o block s of
soil alt hough t he details of t he g r ain mo tion near t he boundar y betw een t he
tw o slices will be v er y much more com plicated wit hout a glue matrix t o hold
t hem in place.
What is t he acceler ation of t he upper slice? Consider t he firs t ε ≪ 1
units of time. Due t o t he v er y s tiff nature of individual g r ains w e will ignore
t he shear def or mation of t he g r ains t hemsel v es. When t = O( ε) t he contact
f orces betw een t he lo w er slice and t he upper slice are able t o tr ansf er suffi-
cient momentum t hat t he upper slice acceler ates lik e a
u
( t) = t t og et her wit h
t he lo w er slice. Ev entuall y as t g ets bigg er so t hat ε ≪ t < 1 t he g r ains at
t he interf ace betw een t he tw o slices begin t o slip, tur n, and rearr ang e as t he
contact f orces ex ceed t he a v ailable contact s trengt h. This occurs no t all at
once but r at her each of t he man y contacts break s at some r andom time so
t hat t he acceler ation of t he upper slice is s till relativ el y smoo t h on t he O( 1)
scale. V ie w ed under t he microscope of v er y short time scales and v er y small
acceler ation de viations, t he acceler ation cur v e is actuall y rough, so t hat w e
90
7.2. A sliding t hought experiment
1
acceleration
time = alower
Acceleration of slices
vs time
(various frictional pro files)
Upper Slice:
di fferent frictional
properties
ε
1
F igure 7.2: A cceler ation of sliding block s of soil in unidirectional loading:
The lo w er slice is cons tr ained t o acceler ate linear l y in time, t he upper slice
reaches an asym p t o tic sliding regime at acceler ation 1. Inter mediatel y , t he
upper slice ma y approach t he asym p t o te from abo v e or belo w , depending
on t he configur ation. F or small times O( ε) where acceler ation is small, t he
block s are lock ed. During liquef action, t he upper slice mus t be loaded be-
y ond t he lock ed regime t o at leas t initiate sliding and rearr ag ement so t hat
φ can chang e.
91
7.2. A sliding t hought experiment
can model it as an s-continuous function.
As time t = O( 1) t he demand placed on t he contacts incre ases t o t he point
where more and more of t hem begin t o slip and as t he y slip t he y tr ansf er
f orces t o o t her s trong er contacts, some of which will f ail due t o t he higher
demand. Since t he acceler ation of t he upper slice is less t han t he acceler ation
of t he lo w er slice, as time goes on t he v elocity of t he upper slice lags t he
bo tt om slice more and more.
Ev entuall y , when t≫ 1 t he relativ e v elocity v
l
v
u
is so lar g e t hat in a mi-
croscopic increment of time t
y
dt t he bo tt om slice mo v es relativ e t o t he upper
slice b y a significant fr action of a g r ain diameter ( v
l
v
u
) t
y
dt D . A t t his
point t he momentum tr ansf er is quite dynamic so t hat g r ains are bouncing
and rolling o v er each o t her .
N o w consider ins tead t he lo w er slice acceler ated b y a sinusoidal w a v e
f or m a
l
(
^
t) = δsin( 2π f t) here acceler ation is measured on t he same scale,
and δ < 1 so t hat t he maximum demand on t he inter -slice momentum tr ans-
f er is less t han t he asym p t o tic sliding amount imagined abo v e but po tentiall y
much g reater t han t he demand f or tim e ε . Also, w e assume t hat w e do no t en-
ter t he dynamic bouncing g r ain regime described abo v e. This is t he regime
w e expect t o oper ate in during an eart hquak e at leas t initiall y prior t o t he
f or mation of liquefied soils where t he inter -g r ain f orces are small relativ e t o
t he t o tal s tress, and t he displacements can g ro w v er y lar g e.
The goal of a continuum model f or g r ain mo tion rele v ant t o liquef action
research t hen is t o describe t he tr ansition from a lock ed pseudo-s tatic as-
sembl y t o a semi-mobile assembl y and t he consequences f or t hat e v olution
on t he pore pressure, s trengt h, and mobility of t he soil.
The tr aditional “backbone cur v e” model of s tress v s s tr ain in soil is based
on prescribing r ules t hat appro ximate t he time dependence of t he s tress-
s tr ain beha vior of t he soil t hrough ultimatel y s tatis tical means. An exam ple
of t he com parison of v arious numerical dynamic models t o real data and
t o each o t her is giv en in (S te w art et al. 2008 ) and t hese types of r ules w ere
applied extensiv el y in comaprison t o centrifug e data b y Liang ( 1995 ).
Ho w e v er , if w e wish t o ha v e a micromechanical explanation f or t he be-
ha vior of soil, it is insufficient t o sim pl y prescribe a set of r ules which tak en
as a whole predict t he bulk scale obser v ations, w e w ould lik e also t o ha v e
some sort of causal connection betw een t he micromechanics and t he bulk scale
r ules.
92
7.3. Basic Mechanical Models Com pared
7.3 Basic Mechanical Models Com pared
The pore pressure equation deriv ed in chap ter 4 sho w s t hat in or der f or pore
pressure t o build up, t he porosity mus t decrease in time, so t hat t he w ater
is squeezed b y t he g r ains. In t his section w e will consider dimensional anal-
yses of v arious po tential ph ysical processes which could contribute t o bulk -
scale w a v e-ener gy dissipation, and com pare t hem t o t he scale of dissipation
which is reasonable in actual pr actice, t hereb y hoping t o eliminate sim ple
ph ysical models t hat predict t he wrong scale of ener gy dissipation, and f o-
cus on models t hat produce reasonable dissipation r ate scales which might
be furt her refined in future research.
The Sliding Glass Plate Model
Consider a sinusoidal w a v e of f = 1 Hz tr a v eling t hrough a loose sand wit h
a v er ag e w a v e speed of c = 200 m/s. Consider a lar g e am plitude of maxi-
mum displacement of A
max
= 10 cm. Then t he s tr ain is t he g r adient of t he
displacement:
γ =
d
dx
A
max
sin( 2π f( x= c t))
So t hat t he peak am plitude of s tr ain is 2π f A
max
= c 0: 003 . T w o la y ers o f
g r ains considered as homog enous sheets of material w ould mo v e relativ e t o
each o t her on a v er ag e on t he or der of 0: 003 g r ain diameters. What is t he fr ac-
tional kinetic ener gy dissipation of such sheets wit h a representativ e kinetic
friction coefficient?
W e calculate t he r atio E
dissip
= E
K
f or tw o sheets of glass 5mm t hick slid-
ing o v er each o t her during one cy cle of a 1Hz w a v e under a com pression of
P
atm
= 1 bar = 100 kP a, assuming continuous sliding at a kinetic friction
coefficient of μ = 0: 65 t o arriv e at an a v er ag e ener gy dissipation per cy cle:
E
dissip
=
∫
cycle
μ A V
rel
P
atm
dt
F or our pre vious sinusoidal w a v e and a displacement betw een plates of
5mm w e arriv e at a t o tal of 3: 8J=m
2
. The a v er ag e kinetic ener gy o v er a cy cle
is:
E
K
=
1
t
cycle
∫
cycle
ρ A V
2
= 2dt
93
7.3. Basic Mechanical Models Com pared
F or an a v er ag e of 0: 81J=m
2
, and a r atio of E
dissip
= E
K
4: 6 meaning t hat
w e expect 4.6 times as much ener gy t o be dissipated as is present as kinetic
ener gy in t he w a v e. Clear l y a w a v e can no t propag ate under such enor mous
dissipation, and in reality tw o actual glass plates w ould ins tead be lock ed
t og et her under such vibr ations, and ha v e a much f as ter w a v e-speed.
F ractionall y F ractured Plate Model
Alt hough it is no t an entirel y accur ate descrip tion, an alter nativ e descrip tion
which has appealing aspects of realism is t hat w e ma y treat t he la y ers of
g r ains as if t he y w ere tw o plates which are sliding essentiall y lock ed, ex cep t
at a r andom number of locations, where t he sliding v elocity is high, so t hat
during a cy cle tw o sliding g r ains displace b y appro ximatel y a g r ain diameter .
W e divide t he plates int o a lar g e number of small squares, and imagine
each one is eit her lock ed, or sliding. The number of sliding g r ains is bino-
miall y dis tributed wit h some p . The p v alue is small so t hat t he net a v er ag e
displacement is of t he appropriate or der e v en t hough w e assume during a
cy cle t hat t he displacement of a sliding fr agment is O( D) t he g r ain diame-
ter . W e require t hat t he net o v er all displacement be accor ding t o t he s tr ain
as bef ore γ 0: 003 so t hat t he displacement of tw o slipping plates one di-
ameter apart is 1 diameter , but non-slipping plates displace 0. Theref ore,
p = γ = 0: 003 . The problem here is t hat alt hough w e’ v e replaced a lar g e
fr action of t he plate slipping at small v elocity , b y a small fr action of t he plate
slipping at lar g e v elocity , t hen if v is t he o v er all relativ e v elocity of t he tw o
plates, v= p is t he relativ e v elocity of tw o slipping parts, and onl y p fr action
of t he plate is slipping, t heref ore t he dissipation r ate is p( v= p) F( v= p) , where
F( v= p) is t he frictional f orce betw een tw o sliding fr agments. Since t he p can-
cels in t he coefficie nt, f or t his t o be dependent on p w e mus t ha v e v elocity
dependent friction , and furt her more, t o achie v e a much lo w er , sus tainable r ate
t han bef ore, w e mus t ha v e F( v= p) ≪ F( v) . A t high speeds, friction mus t
decrease subs tantiall y , b y or ders of magnitude. Mos t frictional surf aces do
no t ha v e such dr as tic v elocity dependent decrease so t hat t his model is un-
appealing on ph ysical g rounds.
R olling Wheel Analogy
Ins tead of considering tw o adjacent slices as sliding o v er each o t her in t he
manner of rigid plates, or using a model of fr agmented plates t o cap ture t he
94
7.3. Basic Mechanical Models Com pared
heterog eneity of inter actions, w e ma y consider se v er al g r ains in a chain t hat
crosses t he slice barrier as f or ming a sort of wheel on which t he upper slice
rolls o v er t he lo w er slice. Clear l y t he friction experienced b y a rigid block
rolling on tin y bearing balls o v er t he floor is usuall y much reduced from
t he friction experienced when w e slide t he same block wit hout t he bearing
balls. P er haps such a model can giv e us a better or der -of-magnitude es timate
of dissipation and is t heref ore rele v ant f or future s tudy and refinement.
Consider a wheel f or med b y a g roup of g r ains whose diameter is n D a
small multiple of t he g r ain diameter . Let t he or der of magnitude of t he rela-
tiv e v elocity of t he tw o ends of t he wheel be γ n D f . If tw o wheels are adjacent,
t hen at t heir region of contact, t he y are r ubbing and t he relativ e v elocity is
O( 2 γ n D f) . W e can imagine t he contact area as of or der D
2
, related t o a sin-
gle g r ain size. In a slice of square cross sectional area A t he t o tal number of
slipping locations is related t o
A
( n D)
2
and t heref ore t he t o tal slipping area is
proportional t o D
2 A
( n D)
2
. The dissipativ e po w er per unit area is t heref ore of
or der :
1
A
A D
2
( n D)
2
μ( 2 γ n D f) σ
v
/
1
n
2
F or wheels in t he r ang e of 5 t o 10 g r ains in diameter , if μ is no t s trongl y
v elocity dependent, w e can g et dissipation r ates in t he r ang e of 1= 25 t o 1= 100
of t he amount in t he sliding plate exam ple, t hereb y reducing t he dissipation
r ate per cy cle from about 4 times t he mean kinetic ener gy t o per haps 5-10%
of t he mean kinetic ener gy in a cy cle.
Alt hough a rolling-wheel model is per haps t oo sim ple, t his anal ysis does
sugg es t t hat relativ el y rigid inter actions on lengt hs of around 5 t o 10 g r ain
diameters might cause a reduction in dissipation b y reducing t he sliding
surf ace area as a fr action of t he cross sectional slice area while retaining
a similar friction coefficient t o t hat in v ol v ed in pure sliding. F uture direc-
tions should concentr ate on t he role of t hese meso-scale semi-rigid inter ac-
tions which ha v e been obser v ed in R echenmacher et al. ( 2011 ). Ot her models
in v ol ving t hese meso-scale s tr uctures might include a bending beam anal-
ogy f or f orce chains, wit h frictional dissipation during bending, or a w a v e-
tr apping-scattering model in which w a v e ener gy injected int o tight clus ters
of g r ains bounces around and is dela y ed in propag ating, and at each bounce
t he kinetic ener gy is spread o v er lar g er and lar g er masses while conser ving
momentum, t hereb y dissipating bulk scale ener gy in a manner similar t o t he
95
7.4. The R ole of Entrop y -lik e Quantities in Modeling Sand
model in chap ter 5 .
7.4 The R ole of Entrop y -lik e Quantities in
Modeling Sand
Se v er al aut hors ha v e introduced an entrop y -lik e quantity in models of g r an-
ular materials. F or exam ple, t he v er y succesful application of inf or mation
entrop y t o predict t he dis tribution of contact f orces in r andom g r anular pack -
ings of frictionless spheres b y Ng an ( 2004 ) which w as v erified via experi-
mental measurements b y Chan and Ng an ( 2005 ). R eal sands are com posed
of non-spherical and decidedl y frictional particles, and so direct application
of t hese ideas t o mechanical models of sand is difficult. Ho w e v er , it is useful
t o consider t he connection of bulk -scale obser v ables t o non-unique g r anular
micro-s tates q specified as a position, orientation, tem per ature, and v eloc-
ity of specific g r ains ha ving specific shapes. A t a particular point in time,
perturbations t o a particular arr ang ement of g r ains q ! q + dq will tend
t o produce a ne w arr ang ement. W e ma y consider certain s tatis tics of t hese
arr ang ements, such as φ( q) , t he porosity . There is an enor mous set of pos-
sible s tates consis tent wit h bulk -scale measurements such as porosity , t o-
tal v olume, and bulk s tresses or s tiffnesses. Giv en a s tate q in t he possible
set, and a particular w a v e-perturbation wit h ener gy dE , w e ha v e a tr ansition
q ! q + dq which will im pl y φ( q) ! φ( q + dq) . If w e giv e all q v alues
consis tent wit h t he initial bulk -scale s tatis tical obser v ations equal probabil-
ity , w e can use t he s tatis tical entrop y of t he dis tribution S = ⟨log( p( q))⟩
t o quantify t he size of t he set of possible s tates. If w e consider a quantity of
mechanical ener gy dE put int o t he sys tem as activ ation ener gy f or rearr ang e-
ment, t hen t here is a set q
′
of ne w s tates reachable b y some pat h where t he
t o tal ener gy is E
0
+ dE t hroughout t he pat h, and t his is partitioned int o g r ain
kinetic ener gy , g r a vitational and contact po tential ener gy , surf ace ener gy (if
g r ains are breaking), and t her mal ener gy . If w e imagine e qual probability of
all reachable pat hs, t his will induce a probability dis tribution on t he tr ansi-
tion φ( q)! φ( q + dq) . Alt hough it is quite difficult t o imagine g etting t his
com plicated expression from firs t principles, w e ne v ert heless expect t hat cer -
tain tr ansitions will be more lik el y t han o t hers: t he φ v alue will tend t o w ar ds
t he mos t common φ reached b y all t he p ossible pat hs.
A sim ple analogy is t he g ame of “Jeng a” in which s tandar dized w ooden
96
7.4. The R ole of Entrop y -lik e Quantities in Modeling Sand
block s are s tack ed on t op of each o t her . The initial s tacking is alter nating
la y ers of t hree block s. Allo w able tr ansitions a re t hose where one block is
tak en from an ywhere in t he middle and placed on t op of t he s tack. If t here
are N identical block s, t here are N! s tates in which t he height of t he s tack
is N units (corresponding t o all t he block s s tack ed on t op of each o t her , but
t he block s are identical so labels can be rearr ang ed, w e ignore t he horizontal
ro tations of t he block s). There are also N! initial s tates, where t he block s are
s tack ed 3 per la y er and ha v e height N= 3 . Ho w e v er , as soon as w e mak e our
firs t mo v e, t he height will be N= 3 + 1 and t he number of possible s tates is
N N! corresponding t o t he f act t hat w e can choose an y of t he N block s t o
mo v e t o t he t op. As t he g ame prog resses t he height v aries from N= 3 t o N and
t he number of possible s tates f or a giv en height v aries from a minimum of
N! t o some maximum v alue corresponding t o t he “mos t probable” height of
a unif or ml y r andoml y chosen leg al Jeng a s tate. Considering height as v er y
similar t o porosity in a sand, w e t heref ore expect t hat t here are inter mediate
porosities which are more probable t han eit her more dense or more loose
arr ang ements.
Considering onl y bulk -scale obser v ations o v er a lar g e number of g r ains,
w e w ould expect t hat certain tr a ject ories of φ w ould tend t o occur more of-
ten t han o t hers. Since o v er all g r ain rearr ang ements are small during t he ini-
tial s tag es of liquef action w e ma y be able t o explain t he linear relationship
betw een dissipated ener gy density and pore w ater pressure de v elopment
obser v ed b y Da vis and Berrill ( 2001 ) b y a sim ple T a y lor series ar gument.
Suppose t hat f or small perturbations from a mos t-probable v alue of φ
f or a giv en bulk scale s tress s tate under t he pre v ailing conditions of s tatic
equilibrium, t hat w e can define a “free ener gy” E T S( φ) f or some “me-
chanical tem per ature” T , t hen f or a certain exter nal ener gy in put, w e ma y
reach s tates corresponding t o a certain region of t he free ener gy po tential.
F or small perturbations from φ
′
t he highes t entrop y v alue f or φ w e ma y
ha v e t he possibility t o represent S as a second or der T a y lor series about φ
′
,
S = S
0
S
2
2
( φ φ
′
)
2
. W e can t hen define a g ener alized “entropic res t oring
f orce” which acts on t he φ coor dinate as t he neg ativ e deriv ativ e of t his “free
ener gy” wit h respect t o chang es in t he φ coor dinate, F( φ) = S
2
( φ φ
′
) .
Our vibr at or y perturbations w ould t hen act b y a combination of t his entropic
res t oring f orce, and some r andom f orcing related t o t he details of t he vibr a-
tion. dφ = S
2
( φ φ
′
)+ σ( dE)◦ dQ f or some s t ochas tic dQ and a scale σ( dE)
t hat depends on t he ener gy in put in each increment. P articular l y if φ is ob-
ser v ed on a bulk scale in v ol ving man y sand g r ains, t he scale of t he s t ochas tic
97
7.4. The R ole of Entrop y -lik e Quantities in Modeling Sand
perturbation σ( dE) might be negligible so t hat φ might be expected t o trend
s trongl y t o w ar ds its maximum entrop y v alue and t hen under go onl y v er y
small oscillations around t hat v alue once it is achie v ed. This model w ould
be one w a y t o unders tand t he pre viousl y w ell-described tendency f or soils
t o come t o equilibrium at a “critical s tate” wit h eit her cy clic or continuous
direct shearing of appro ximatel y cons tant am plitude and frequency or r ate.
Alt hough t his model is p er haps t oo sim ple, it could in a s tr aightf or w ar d w a y
be com pared t o experimental data wit h t he appropriate appar atus, and mod-
ified appropriatel y . F uture w or k should k eep in mind t his g r anular micro-
s tate vie wpoint and consider whet her it can be expanded b y considering
o t her s tate v ariables besides φ . F urt her more because during liquef action t he
chang es in φ required t o induce pore pressure are actuall y quite small, t he
model need onl y w or k o v er a small r ang e of φ v alues relativ e t o t he initial
v alue, so t hat t he T a y lor series approach becomes a v alid means of sim plify -
ing t he dynamics.
98
Chapter 8
Conclusion
Desdemona: Oh mos t lame and im po t ent conclusion!
— W illiam Shak espeare Ot hello
By appl ying mat hematical modeling techniques suitable f or g ener al mechan-
ical problems t o specific problems in soil lique f action and t her mal w a v e dis-
sipation t his w or k highlights tw o im portant fundamental ph ysical phenom-
ena: t he flo w of w ater t hrough soil during liquef action which contr adicts t he
undr ained assum p tions of s tandar d liquef action research, and t he dissipa-
tion of w a v e ener gy t hrough t her mall y driv en momentum diffusion which
po tentiall y expla ins t he detectability of v er y tin y fr actures in cr ys talline ma-
terials. These disco v eries pro vide a whole ne w direction f or t he s tudy of
soil liquef action and a s tarting point f or furt her s tudy of nano-scale w a v e
propag ation.
F uture directions f or research include t he possibility of appl ying t he NS A
based continuum modeling approach t o t he deriv ation of continuum mod-
els f or g r ain mo tion. These could couple t o t he fluid flo w equations f or t he
prediction of liquef action in realis tic heterog eneous deposits. In addition,
t he Ba y esian com putational techniques used t o identify t he par ameters of
t he w a v e dissipation model could be utilized t o identify t he par ameters of
g r ain mo tion equations based on obser v ations in eit her real eart hquak es via
boreholes, or in g eo technical centrifug e models.
F urt her in v es tig ation of micro-scale w a v e propag ation and dissipation
could include obser ving t he dis tribution of small magnitude e v ents in table-
t op experiments such as t he s tress fr acturing of a cr ys tal in Ås tröm et al.
( 2006 ) and pseudo-s tatic sliding on a simulated f ault in V oisin, R enar d, and
Gr asso ( 2007 ). By specifying a t hreshold f or detection of a w a v e, w e can es ti-
99
mate t he effect of small-signal censoring and deter mine ho w t his censoring
of v er y small e v ents affected t he obser v ed his t og r am of e v ent sizes.
In t he context of Ba y esian modeling of dynamical sys tems, t he difficul-
ties t hat arise in doing efficient inf erence whene v er com putational models
are used t o obtain predictions sugg es ts an area of research in v ol ving making
model errors explicit t hrough a choice of k er nel or pseudo-k er nel such as t he
Gaussian process model utilized here. In t he presence of such a modeling-
error model per haps t he continuity of t he models can be used t og et her wit h a
fix ed database of r uns t o inter polate t he errors f or inter mediate v alues of im-
portant par ameters, and obtain better inf erence wit hout enor mous amounts
of com putation. F uture research could include sof tw are packag es t hat aut o-
maticall y use t hese techniques t o pro vide f as t efficient inf erence f or dynamic
models.
100
Bibliograph y
Ås tröm, J., P .C.F . Di S tef ano, F . Pröbs t, L. S t odolsky , J. T imonen, C. Bucci,
S. Cooper , et al. 2006. “F r acture processes obser v ed wit h a cr y og enic
detect or . ” Phy sics Lett er s A 356, no. 4–5 (A ugus t): 262–266. :
0375-9601, accessed Januar y 25, 2013.
doi: 10 . 1016 /j.physleta. 2006 . 03 . 059 . http://www.sciencedirect.
com/science/article/pii/S 037 5960106005007 .
Barenblatt, G. I. 2003. Scaling. Cambridg e texts in applied mat h ematics.
Cambridg e: Cambridg e U niv ersity Press. : 0521826578.
Bear , Jacob. 1988. Dynamics of fluids in por ous media. Do v er book s on ph ysics
and chemis tr y . N e w Y or k: Do v er . : 0486656756.
Bear d, D.C., and P .K. W e y l. 1973. “Influence of T exture on P orosity and
P er meability of U nconsolidated Sand. ” AAPG Bulletin 57. :
0149-1423, accessed Sep tember 18, 2011.
doi: 10 . 1306 / 819 A 427 2 - 16 C 5 - 11 D 7 - 8645000102 C 1865 D .
http://search.datapages.com/data/doi/ 10 . 1306 / 819 A 427 2 - 16 C 5 -
11 D 7 - 8645000102 C 1865 D .
Brinkman, H. C. 1949. “ A calculation of t he viscous f orce ex erted b y a
flo wing fluid on a dense sw ar m of particles. ” Applied Scientific Resear c h 1
(December): 27–34. : 0365-7132, 0003-6994, accessed Sep tember 19,
2011. doi: 10 . 1007 /BF 02120313 .
http://www.springerlink.com/index/ 10 . 1007 /BF 02120313 .
Chan, S. H., and A . H. W . Ng an. 2005. “S tatis tical dis tribution of contact
f orces in packings of def or mable spheres. ” Mec hanics of Mat erials 37 (4):
493–506. : 0167-6636.
101
Bibliog r aph y
Da vis, R. O., and J. B. Berrill. 2001. “P ore Pressure and Dissipated Ener gy in
Eart hquak es—F ield V erification. ” Journal of Geo t ec hnical and
Geoen vir onment al Engineering 127 (3): 269. : 10900241.
doi: 10 . 1061 /(ASCE) 1090 - 0241 ( 2001 ) 127 : 3 ( 269 ) .
http://link.aip.org/link/JGGEFK/v 127 /i 3 /p 269 /s 1 &Agg=doi .
F ieg el, Gregg L., and Br uce L. K utter . 1994. “Liquef action Mechanism f or
La y ered Soils. ” Journal of Geo t ec hnical Engineering 120 (4): 737. :
07339410. doi: 10 . 1061 /(ASCE) 07 33 - 9410 ( 1994 ) 120 : 4 ( 7 37 ) .
http://link.aip.org/link/JGENDZ/v 120 /i 4 /p 7 37 /s 1 &Agg=doi .
F igueroa, J. L udwig, A del. S. Saada, Liqun Liang, and Nitin M. Dahisaria.
1994. “Ev aluation of Soil Liquef action b y Ener gy Principles. ” Journal of
Geo t ec hnical Engineering 120 (9): 1554–1569.
http://dx.doi.org/ 10 . 1061 /(ASCE) 07 33 - 9410 ( 1994 ) 120 : 9 ( 1554 ) .
F o w ler , Andre w . 1998. Mat hematical models in t he applied sciences. (R epr .)
Cambridg e: Cambridg e U niv . Press. : 9780521467032.
Gal v anett o, Ugo, and M. H. Aliabadi, eds. 2010. Multiscale modeling in solid
mec hanics: com put ational appr oac hes. Com putational and experimental
met hods in s tr uctures, v ol. 3. London : London ; N e w Y or k: Im perial
Colleg e ; Dis tributed b y W or ld Scientific. : 9781848163072.
Ge y er , Char les J., and Leif Johnson. 2013. mcmc: Mar kov Chain Mont e Car lo,
March. A ccessed Jul y 30, 2013.
http://cran.r- project.org/web/packages/mcmc/index.html .
Goren, L., E. Aharono v , D. Spar k s, and R. T oussaint. 2010. “P ore pressure
e v olution in def or ming g r anular material: A g ener al f or mulation and
t he infinitel y s tiff appro ximation. ” J. Geophy s. Res 115.
. 2011. “ The Mechanical Coupling of Fluid-F illed Gr anular Material
U nder Shear . ” Pur e and Applied Geophy sics (June). : 0033-4553,
accessed Sep tember 9, 2011. doi: 10 . 1007 /s 00024 - 011 - 0320 - 4 .
http://www.springerlink.com/index/ 10 . 1007 /s 00024 - 011 - 0320 - 4 .
102
Bibliog r aph y
Higdon, Da v e, Marc K ennedy , James C. Ca v endish, John A . Caf eo, and
R obert D. R yne. 2004. “Combining F ield Data and Com puter
Simulations f or Calibr ation and Prediction. ” SIAM Journal on Scientific
Com puting 26, no. 2 (Januar y): 448–466. : 1064-8275, 1095-7197,
accessed Jul y 29, 2013. doi: 10 . 1137 /S 1064827 503426693 . http://epubs.
siam.org.libproxy.usc.edu/doi/abs/ 10 . 1137 /S 1064827 503426693 .
Holtz, R, and W K o v acs. 1981. An intr oduction t o g eo t ec hnical engineering.
Engle w ood Cliffs N .J.: Prentice-Hall. : 9780133539127.
Ho wison, Sam. 2005. Pr actical applied mat hematics: modelling, anal y sis,
appr o ximation. Cambridg e texts in applied mat hematics. N e w Y or k:
Cambridg e U niv ersity Press. : 0521842743.
Idriss, I. M., and R. W . Boulang er . 2006. “Semi-em pirical procedures f or
e v aluating liquef action po tential during eart hquak es. ” Soil Dynamics and
Eart hq uake Engineering 26 (2-4): 115–130.
Kilmins ter , De vin, and R eason Machete. 2005. “Prediction, Beha viour , and
Ignor ance. ” In Pr oceedings of 2005 Int ernational Sym posium on N onlinear
Theor y and its Applications. Br ug es, Belgium. : 4-88552-215-3, accessed
Jul y 26, 2013. http:
//www.ieice.org/proceedings/NOLTA 2005 /HTMLS/paper/ 7 115 .pdf .
K okusho, T ak eji. 1999. “W ater F ilm in Liquefied Sand and Its Effect on
Later al Spread. ” Journal of Geo t ec hnical and Geoen vir onment al Engineering
125 (10): 817. : 10900241.
doi: 10 . 1061 /(ASCE) 1090 - 0241 ( 1999 ) 125 : 10 ( 817 ) .
http://link.aip.org/link/JGGEFK/v 125 /i 10 /p 817 /s 1 &Agg=doi .
K okusho, T ak eji, and T etsuro K ojima. 2002. “Mechanism f or
P os tliquef action W ater F ilm Gener ation in La y ered Sand. ” Journal of
Geo t ec hnical and Geoen vir onment al Engineering 128, no. 2 (F ebr uar y):
129–137. http://link.aip.org/link/?QGT/ 128 / 129 / 1 .
Kr amer , S te v en. 1996. Geo t ec hnical eart hq uake engineering. U pper Saddle Riv er
N .J.: Prentice Hall. : 9780133749434.
Lak eland, Daniel L., Am y R echenmacher , and R og er Ghanem. n.d.
“ T o w ar ds a Com plete Model of Soil Liquef action: t he im portance of
fluid flo w and g r ain mo tion. ” Pr oceedings of t he Roy al Society of London.
Series A , Mat hematical and Phy sical Sciences.
103
Bibliog r aph y
Liang, Liqun. 1995. “De v elopment of an ener gy met hod f or e v aluating t he
liquef action po tential of a soil deposit. ” PhD diss., Case W es ter n
R eser v e U niv ersity .
http://rave.ohiolink.edu/etdc/view?acc_num=case 1058541489 .
Lobr y , C., and T . Sari. 2008. “N on-s tandar d anal ysis and representation of
reality . ” Int ernational Journal of Contr ol 81 (3): 519–536. : 0020-7179.
http://www.informaworld.com/ 10 . 1080 / 00207 17 07 016017 28 .
MacKa y , Da vid. 2003. Inf ormation t heor y , inf er ence, and learning algorit hms.
Cambridg e UK ;;N e w Y or k: Cambridg e U niv ersity Press. :
9780521642989.
Maha jan, Sanjo y . 2010. S tr eet-fighting mat hematics: t he art of educat ed guessing
and opportunis tic pr oblem sol ving. Cambridg e, Mass: MIT Press. :
9780262514293.
Ma jmudar , T . S., and R. P . Behring er . 2005. “Contact f orce measurements
and s tress-induced aniso trop y in g r anular materials. ” N atur e 435, no.
1079 (June): 1079–1082. : 0028-0836. doi: 10 . 1038 /nature 03805 .
http://www.nature.com/doifinder/ 10 . 1038 /nature 03805 .
Marin, Jean-Michel, Pierre Pudlo, Chris tian P . R obert, and R obin J. R yder .
2012. “ Appro ximate Ba y esian com putational met hods” [in en]. S t atis tics
and Com puting 22, no. 6 (N o v ember): 1167–1180. : 0960-3174,
1573-1375, accessed Jul y 26, 2013. doi: 10 . 1007 /s 11222 - 011 - 9288 - 2 .
http:
//link.springer.com.libproxy.usc.edu/article/ 10 . 1007 /s 11222 -
011 - 9288 - 2 .
N elson, E. 1977. “Inter nal set t heor y : a ne w approach t o nons tandar d
anal ysis. ” AMERICAN MA THEMA TICAL SOCIET Y 83 (6): 1165.
N emat-N asser , S., and A . Shok ooh. 1979. “ A unified approach t o
densification and liquef action of cohesionless sand in cy clic shearing. ”
Canadian Geo t ec hnical Journal 16 (4): 659–678. : 1208-6010.
doi: 10 . 1139 /t 7 9 - 07 6 . http://www.nrc.ca/cgi-
bin/cisti/journals/rp/rp 2 _abst_e?cgj_t 7 9 -
07 6 _ 16 _ns_nf_cgj 16 - 7 9 .
104
Bibliog r aph y
Ng an, A . H. W . 2004. “On dis tribution of contact f orces in r andom g r anular
packings. ” Phy sica A: S t atis tical Mec hanics and its Applications 339 (3-4):
207–227. : 0378-4371.
NIS T . 2012. Thermophy sical Pr operties of Fluid Sy s t ems. A ccessed F ebr uar y 14.
http://webbook.nist.gov/chemistry/fluid/ .
Nyquis t, Chris tina. 2012. The Mar c h 11 T ohok u Eart hq uake, One Y ear Lat er . What
Ha v e W e Learned? | Science F eatur es. A ccessed Januar y 11, 2013. http:
//www.usgs.gov/blogs/features/usgs_top_story/the- march- 11 -
tohoku- earthquake- one- year- later- what- have- we- learned/ .
Plim p t on, S te v e. 1995. “F as t P ar allel Algorit hms f or Short-Rang e Molecular
Dynamics. ” Journal of Com put ational Phy sics 117, no. 1 (March): 1–19. :
0021-9991, accessed F ebr uar y 1, 2013. doi: 10 . 1006 /jcph. 1995 . 1039 .
http://www.sciencedirect.com/science/article/pii/
S 00219991857 1039 X .
R echenmacher , Am y L., Sar a A bedi, Olivier Chupin, and
Andrés D. Or lando. 2011. “Char acterization of mesoscale ins tabilities in
localized g r anular shear using digital imag e correlation” [in en]. A ct a
Geo t ec hnica 6, no. 4 (December): 205–217. : 1861-1125, 1861-1133,
accessed A ugus t 19, 2013. doi: 10 . 1007 /s 11440 - 011 - 0147 - 2 . http:
//link.springer.com.libproxy.usc.edu/article/ 10 . 1007 /s 11440 -
011 - 0147 - 2 .
R obert, Alain. 2011. N ons t andar d anal y sis. Mineola N .Y .: Do v er Publications.
: 9780486432793.
R obinson, A br aham. 1996. N on-s t andar d anal y sis. R e v . ed. Princet on N .J.:
Princet on U niv ersity Press. : 9780691044903.
Sa wicki, Andr zej, and Jecek Mierczyński. 2006. “De v elopments in Modeling
Liquef action of Gr anular Soils, Caused b y Cy clic Loads. ” Applied
Mec hanics Review s 59 (2): 91. : 00036900. doi: 10 . 1115 / 1 . 2130362 .
http://link.aip.org/link/AMREAD/v 59 /i 2 /p 91 /s 1 &Agg=doi .
Seed, R.B., K.O. Cetin, R.E.S. Moss, A .M. Kammerer , J. W u, J.M. P es tana,
M.F . Riemer , R.B. Sancio, J.D. Br a y , R.E. Ka y en, et al. 2003. “R ecent
adv ances in soil liquef action engineering: a unified and consis tent
fr ame w or k. ” In Pr oceedings of t he 26t h Annual ASCE Los Ang eles
Geo t ec hnical Spring Seminar : Long Beac h, CA.
105
Bibliog r aph y
Siv at ha y alan, S ., and Y .P . V aid. 1998. “ T r ul y undr ained response of g r anular
soils wit h no membr ane-penetr ation effects. ” Canadian Geo t ec hnical
Journal 35 (5): 730–739. : 12086010. doi: 10 . 1139 /cgj- 35 - 5 - 7 30 . http:
//www.nrc.ca/cgi- bin/cisti/journals/rp/rp 2 _abst_e?cgj_t 98 -
048 _ 35 _ns_nf_cgj 35 - 98 .
Snelson, Edw ar d, Car l Edw ar d Rasmussen, and Zoubin Ghahr amani. 2004.
“W ar ped g aussian processes. ” Adv ances in neur al inf ormation pr ocessing
sy s t ems 16:337–344. A ccessed A ugus t 2, 2013.
http://books.google.com/books?hl=en&lr=&id= 0 F- 9 C 7 K 8 fQ 8 C&oi=
fnd&pg=PA 337 &dq=% 22 a+form+well+modelled+by+a+GP.+Making+
such+a+transformation+should+really+be+a% 22 +% 22 and+
asymmetric+noise+in+general.+It+is+not+however+just+a+GP+
with+a% 22 +&ots=TGKvoYRa 50 &sig= 7 CqnWm 3 _OWmbDU 6 iAg 1 cHPzJeZo .
S te w art, J. P ., A . O. L. Kw ok , Y . M. A . Hashash, N . Mataso vic, R. Pyk e,
Z. W ang, and Z. Y ang. 2008. “P acific Eart hquak e Engineering R esearch
Center . ” U niv er sity of Calif ornia, Ber kele y. A ccessed Januar y 22, 2013.
http:// 128 . 32 . 145 . 86 /publications/peer_reports/reports_ 2008 /
web_PEER 804 _STEWARTetal.pdf .
T addy , Matt he w A ., Herbert K. H. Lee, and Br uno Sansó. 2009. “F as t
inf erence f or s tatis tical in v erse problems” [in en]. In v er se Pr oblems 25, no.
8 (A ugus t): 085001. : 0266-5611, accessed Jul y 31, 2013.
doi: 10 . 1088 / 0266 - 5611 / 25 / 8 / 085001 .
http://iopscience.iop.org/ 0266 - 5611 / 25 / 8 / 085001 .
The Maxima Project. 2011. Maxima, A Com put er Alg ebr a Sy s t em. V er sion 5.25.0.
http://maxima.sour cef or g e.net. http://maxima.sourceforge.net .
V erschoor , H. 1951. “Experimental data on t he viscous f orce ex erted b y a
flo wing fluid on a dense sw ar m of particles. ” Applied Scientific Resear c h 2
(Januar y): 155–161. : 0365-7132, 0003-6994, accessed Sep tember 21,
2011. doi: 10 . 1007 /BF 0041197 9 .
http://www.springerlink.com/index/ 10 . 1007 /BF 0041197 9 .
V oisin, Chris t ophe, F r ançois R enar d, and Jean-R obert Gr asso. 2007. “Long
ter m friction: F rom s tick -slip t o s table sliding” [in en]. Geophy sical
Resear c h Lett er s 34 (13): n/a–n/a. : 1944-8007, acce ssed Januar y 25,
2013. doi: 10 . 1029 / 2007 GL 0297 15 . http:
//onlinelibrary.wiley.com/doi/ 10 . 1029 / 2007 GL 0297 15 /abstract .
106
Bibliog r aph y
Zienkie wicz, O. C., C.T . Chang, and P . Bettess. 1980. “Dr ained, undr ained,
consolidating and dynamic beha viour assum p tions in soils. ”
Géo t ec hniq ue 30, no. 4 (Januar y): 385–395. : 0016-8505, 1751-7656,
accessed Sep tember 3, 2013. doi: 10 . 1680 /geot. 1980 . 30 . 4 . 385 .
http://www.icevirtuallibrary.com.libproxy.usc.edu/content/
article/ 10 . 1680 /geot. 1980 . 30 . 4 . 385 .
107
Abstract (if available)
Abstract
The research in this dissertation encompasses four main thrusts: mathematical modeling techniques for building continuum models of systems that are explicitly acknowledged to be non-continuous at the finest scale of observation, application of a continuum model of water flow in porous media to the problem of “drained vs undrained” processes in soil liquefaction, an example model for the dissipation of waves in a system of molecules decoupled from any radiative process, where the total energy stays constant, and the application of Bayesian Statistics to identifying the parameters of these physical models.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Uncertainty management for complex systems of systems: applications to the future smart grid
PDF
Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection
PDF
A stochastic Markov chain model to describe cancer metastasis
PDF
Feature and model based biomedical system characterization of cancer
PDF
Physics-based bistatic radar scattering model for vegetated terrains in support of soil moisture retrieval from signals of opportunity
PDF
Inverse modeling and uncertainty quantification of nonlinear flow in porous media models
PDF
The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems
PDF
Sense and sensibility: statistical techniques for human energy expenditure estimation using kinematic sensors
PDF
Design optimization under uncertainty for rotor blades of horizontal axis wind turbines
PDF
Developing frameworks to quantify the operational and environmental performance of energy systems within the context of climate change
PDF
Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties
PDF
New approaches to satellite formation-keeping and the inverse problem of the calculus of variations
PDF
Cache analysis and techniques for optimizing data movement across the cache hierarchy
PDF
Electrokinetic transport of Cr(VI) and integration with zero-valent iron nanoparticle and microbial fuel cell technologies for aquifer remediation
Asset Metadata
Creator
Lakeland, Daniel L.
(author)
Core Title
Continuum modeling techniques and their application to the physics of soil liquefaction and dissipative vibrations
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Civil Engineering
Publication Date
11/18/2014
Defense Date
10/23/2013
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
Bayesian statistics,continuum models,energy dissipation,nonstandard analysis,OAI-PMH Harvest,soil liquefaction,statistical mechanics,vibrations
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Ghanem, Roger G. (
committee chair
), Ben-Zion, Yehuda (
committee member
), Lynett, Patrick J. (
committee member
), Newton, Paul K. (
committee member
), Rechenmacher, Amy Lynn (
committee member
)
Creator Email
dlakelan@street-artists.org,lakeland@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-346615
Unique identifier
UC11297463
Identifier
etd-LakelandDa-2156.pdf (filename),usctheses-c3-346615 (legacy record id)
Legacy Identifier
etd-LakelandDa-2156.pdf
Dmrecord
346615
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Lakeland, Daniel L.
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
Bayesian statistics
continuum models
energy dissipation
nonstandard analysis
soil liquefaction
statistical mechanics
vibrations