University of Southern California Dissertations and Theses

Ice: an impedance spectroscopy and atomistic simulation study

(USC Thesis Other)

Ice: an impedance spectroscopy and atomistic simulation study

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

Ice: An Impedance Spectroscopy and
Atomistic Simulation Study

A PhD Dissertation Manuscript Submitted by:

Keith Chin

Mork Family Department of Chemical Engineering and
Materials Science
University of Southern California
June 2014

i

Acknowledgements
Obviously, my PhD work is the final culmination of many years of aggravation,
perseverance, and dedication to a fundamental belief that education is the foundation to success
in life. However, I would not have completed this herculean effort without many individuals’
support. The most important people I would like to thank are my family, especially my wife
Victoria. I don’t know how she put up with me through all the years of helplessly watching me
dig into the visceral of my research and overcoming one problem after another. I only hope she
forgives me for those years of neglect as my mind incessantly wonders away during our
conversions thinking only of work. To this, I promise to make it up to her for all those years
lost. To my daughter, Loryn, I hope to gain back all those precious time lost as well. You are an
amazing young lady even at the age of 13. For better or worse, I hope my dissertation work will
also encourage you to pursue all your goals in life despite any level of difficulty. I can already
see that you are by far more intelligent than me but you’ll have to learn to endure and never give
up on your dreams and goals like your father. Just remember, as long as you have your health
you can do anything. As for my parents, I only hope they truly understand what this
accomplishment means to our family of first generation Chinese immigrants. I know it must be
difficult to make ends meet after migrating to a new country to improve our lives. Hopefully,
now you see what truly this great country has to offer in terms of providing equal opportunity for
all who is willing to work hard and to strive for a better life beyond just monetary gains.
My employment at the Jet Propulsion Laboratory also had a large part in my PhD work.
The tuition reimbursement program at JPL was a big help financially and provided me with an
opportunity to compete my PhD work without placing a huge financial burden on my family. To
this, I will strive to perform my very best in supporting all NASA goals for the future and uphold
all the integrity that this great institution represents. Many colleagues at JPL also supported this
endeavor through the years, especially Dr. Martin Beuhler who initiated my dissertation
experimental work at JPL. I hope you are finally enjoying your retirement with your family in
the great Palouse of Washington state. To Dr. Bob Herman and his wife and good friend
Jennifer Herman, I thank you guys for being so supportive through both good and bad times. It’s

ii

always easier to strive for something so challenging with good friends and colleagues behind you
with words of encouragement.
Finally, to all my friends and colleague at USC I really appreciate all your help beyond
words can convey. Dr. Wei Tao is one of the best scientists I know even before he graduated.
His diligence and expertise in his research fields is second to none. I only hope I can come
remotely close to your level of comfort and understanding of simulations someday. To Size
Zheng, I really appreciate all your assistance as well. I have the highest confidence you will
become a great scientist soon. You have some amazing talents already. To my advisor, Dr.
Katherine Shing, your efforts to make me a better scientist may someday prove worthwhile. For
now I just relieved to complete what I’ve set out to do years and years back. I only wish I could
give you the “perfect” dissertation, but given all my time limitations it is a major challenge
obviously. I know you’ve understood my circumstances well and I really appreciate your
patience throughout the process. Despite it all, I will heed your words and will strive to be as
through as possible in all my scientific endeavors.

iii

Table of Contents

Chapter 1: Introduction ----------------------------------------------------------------------------- 1

Chapter 2: Equipment, Experiment, and Calibration ---------------------------------- 6
2.1 Background ------------------------------------------------------------------------------------------ 6
2.2 Experimental----------------------------------------------------------------------------------------- 6
2.3 Calibration Results & Discussion --------------------------------------------------------------- 11
2.3.1 Instrument and Measurement Calibration --------------------------------------------- 11
2.3.2 Calibration on Electrode Configuration with Pure DI Water ---------------------- 14
2.3.3 Calibration on Liquid Electrolyte Measurement Results ---------------------------- 17
2.4 Conclusions ---------------------------------------------------------------------------------------- 24

Chapter 3: Electrical Spectroscopy of Ice ------------------------------------------------- 25
3.1 Theory ---------------------------------------------------------------------------------------------- 25
3.1.1 The Debye Theory for Relaxation in Ice ----------------------------------------------- 27
3.1.2 Polarization and Resonance in Ice from a Molecular Perspective ----------------- 29
3.1.3 Conduction in ice ------------------------------------------------------------------------- 31
3.2 Literature Review on Ice Electrical Properties ------------------------------------------------ 33
3.3 In-Situ Ice Observations -------------------------------------------------------------------------- 37
3.4 Measurement Results on Ice --------------------------------------------------------------------- 41
3.4.1 Ice Doped with KCl ---------------------------------------------------------------------- 43
3.4.2 Ice Doped with KOH --------------------------------------------------------------------- 45
3.4.3 Ice Doped with CaCl
2
-------------------------------------------------------------------- 45
3.4.4 Ice Doped with KI, NaCl, and MgSO
4
-------------------------------------------------- 46
3.5 Equivalent Circuit Modeling on Ice ------------------------------------------------------------- 65
3.5.1 The Debye Equivalent Circuit --------------------------------------------------------- 65
3.5.2 The Debye Equivalent Circuit with Leakage ----------------------------------------- 68
3.5.3 The Debye Equivalent Circuit with Proton Conduction and Resonance --------- 72
3.5.4 Equivalent Circuit Modeling Results -------------------------------------------------- 76
3.6 Results and Discussion -------------------------------------------------------------------------- 84

iv

2.6.1 Debye Relaxation Times ----------------------------------------------------------------- 84
2.6.2 Debye Resonance Times ----------------------------------------------------------------- 90
2.6.3 Ice Conduction ---------------------------------------------------------------------------- 94
3.7 Summary and Conclusion --------------------------------------------------------------------- 102

Chapter 4: Molecular Dynamic Simulations on Ice ---------------------------------- 104
4.1 Literature Review on Ice Simulations -------------------------------------------------------- 104
4.2 Theory and Simulation Details ----------------------------------------------------------------- 110
4.3 Results and Discussion ------------------------------------------------------------------------- 118
4.3.1 Kinetic Energy, Potential Energy, and Density ------------------------------------- 118
4.3.2 Structural Properties from the Radial Distribution Function --------------------- 128
4.3.3 Ionic Transport Mechanism and Jump Events --------------------------------------- 135
4.3.4 Ionic Transport Properties and the Mean-Square Displacement ----------------- 148
4.3.5 Electrical Properties -------------------------------------------------------------------- 155
4.3.6 Visualization and Trajectory Plots ---------------------------------------------------- 162
4.4 Conclusions --------------------------------------------------------------------------------------- 169

Chapter 5: Final Conclusions ------------------------------------------------------------------ 172
5.1 Suggested Future Work ------------------------------------------------------------------------- 176

Appendix A: Fundamental Equations for Conduction & Polarization by AC Impedance ----- 177
Appendix B: Derivation of the Debye Equivalent Circuit from Basic Principles --------------- 179
Appendix C: Equivalent Circuit Model Derivations ------------------------------------------------ 184
C.1: Derivation of the Debye Equivalent Circuit ------------------------------------------------ 184
C.2: Derivation of the Debye Equivalent Circuit with Leakage ------------------------------- 186
C.3: Derivation of the Debye Equivalent Circuit with Leakage and Resonance ------------ 188

Symbols ---------------------------------------------------------------------------------------------------- 191
Bibliography ---------------------------------------------------------------------------------------------- 195
References ------------------------------------------------------------------------------------------------- 196

1

1.0 Introduction
H
2
O has been the subject of one of the most extensively studied substance due to its
importance in nature and all living organisms on Earth. Bulk ice and ice particles also play a
basic role in terrestrial, planetary, and atmospheric sciences on Earth. As mankind approaches
the cusp of deep space exploration in our quest to seek out new or similar forms of life, ice has
experienced a resurgence of new research interests due to its omnipresence beyond Earth.
Though existence of liquid water is yet to be determined, ice is known to exist on virtually all
celestial bodies.
1,2,3,4
Therefore, the search for life beyond our planet hinges on the
understanding of planetary ice and its interactions with bio-precursors that could preserve
indications of pre-existent life or currently capable of sustaining ecosystems suitable for life by
providing nutrients in the form of chemical energy.
To enable life detection in future spacecraft landers and rovers, advanced in-situ
instrumentation development holds the key to providing unambiguous identification of water in
future exploratory mission to icy planets. Such an instrument should be capable of not only
identifying the presence of water, but also determine the phase, compositional constituents, and
quantification of water in both homogeneous and heterogeneous mixtures containing water. Up
to now, most remote sensing flight instruments designed for observation of water is based on
inferences from indirect measurements (e.g., images, spectroscopic reflection, orbital motion,
etc.).
5
For example in arid soils, there is currently no way to directly measure water thin films or
their chemistry. Knowing the presence of film water would address one of the most difficult and
important questions in frozen-arid environments: How do ions (and, for earth, microbes) get
transported within the soils under extreme environmental conditions? Cryoturbation is one
possibility, but thin film transport is also possible, especially if brines are present.
6
The most
viable and unambiguous identification of water necessitates an in-situ analysis of distinctive
chemical, physical, or thermodynamic properties of the individual water molecule or the
aggregate behavior of water molecules in response to an external stimulus relative to the existing
environmental conditions (e.g., thermodynamic, physical, morphological, etc.) conditions. Thus,
measurements involving electromagnetic, thermal, or mechanical techniques offer the most
reliable means of providing direct stimulation and measureable responses for potential
application in life detection in water in both liquid and ice phases.

2

Though several in-situ methods exist or are in development for possible life detection
under extreme environments, ac electrochemical impedance spectroscopy (EIS or IS) offers the
most advantages. Development of an EIS system for space applications is relatively simple, and
therefore, very cost effective with scalability down to IC-level in size. There are also numerous
possibilities of deployment platforms associated with EIS (e.g., underground drills, robotic arms,
insertion probes, etc.). Unlike techniques involving sample heating and analyzing thermally
evolved volatiles (e.g., chromatographic methods)
7
, EIS offers non-destructive measurement and
evaluation of samples under investigation and thereby provide unlimited sampling capabilities, at
least for relatively accessible samples. DC electrical and thermal conductivity measurements are
seriously limited by porosity issues in heterogeneous mixtures containing hydrated minerals and,
thus, render them incapable of distinguishing between water and ice measurements; electrical
spectroscopy techniques offer sensitivity between water and ice while still retaining the
advantages from conductivity probes.
8,9

All dipole molecules exhibit a strong and unique characteristic response from EIS
measurements due to the frequency dependent relaxation and resonance processes emanating
from the molecular level. At high frequencies (MHz), the response is due to atomic polarization
caused by the stretching of individual bonds within the water molecule. Due to nature of the
hydrogen network of ices, the well-known characteristic response is the slow relaxation process
associated with orientational motion. Ice also exhibits the double-layer response similar to liquid
water at the electrode interface due to charge balance considerations. In soils containing water
and/or ice, the characteristic response is very similar to that of pure ice as dictated by the
water/ice content. In these heterogeneous mixture of materials, relaxation processes will only
occur with the existence of bulk water or ice at high water contents; thus, IS responses at low
water/ice content is dominated by the electrical behavior of soil grains since the water and ice
only exist in the form of interfacial films at the grain boundaries.
10
The interpretation of IS
spectra is further complicated by the concentration of charge carriers present in the water solvent
due to ionic transport processes in response to the applied electric field. The dependency on
ionic concentration can also be readily observed by change in resonance and relaxation
frequencies as well as the effective conductivity of the bulk material. Therefore, detailed
understanding of solute-solvent interactions is the most arduous hurdle in the unambiguous
interpretation and characterization of in-situ IS responses.

3

Design and development of an instrument capable of ac impedance analysis also pose
significant challenges in terms of retaining sensitivity requirements along with durability under
field conditions. Historically, almost all experimental data on the dielectric behavior of ice were
conducted using parallel plate electrode configuration back.
11,12
The parallel contact electrodes
will not adequately endure durability requirements in subsurface measurements and often lead to
ice fractures in in-situ measurements under extreme thermodynamic conditions.
13,14
Even non-
uniform point electrodes results in “imperfect” contacts to solid electrolytes which has a
profound influence on the frequency dependence of the potential distribution at the interface.
15

Irrespective of electrode geometrical issues, electrode polarization is also frequency dependent
and constitutes a major consideration for interpretation of EIS measurements as well.
16

Since EIS is a bulk measurement technique, it can only provide an average or an effective
characterization of the sample material and cannot fully probe the complex physical and
chemical mechanisms emanating at the nanoscale level. Modern theories on ice transport
mechanisms is still mostly inferred from the classical thought based on the formation and
migration of protonic defects, more commonly known as ionic and Bjerrum defects. It is
generally understood that ionic defects are independent of temperature and are formed by the
dissociation of water molecules into OH
-
and H
3
O
-
ions by a process known as tunneling where
the proton from adjacent water molecules hop only the bond between them. On the contrary,
Bjerrum defects (both L,D-types) have been the focus of many molecular studies since they are
temperature dependent, and thus, exhibit much more dynamic behavior.
17,18,19
At best, results
from these early works were very vague in concluding that the L,D-pair dynamics is a very
complex process involving local melting of ice and H-bond changes in groups of molecules.
20

Additionally, the emphasis on protonic defects from these past molecular studies are typically
limited due to simulation of relatively small clusters (<100 molecules) which pose substantial
challenges in energy minimization (allowing the displacement of the molecules) techniques
thereby rending them highly inaccurate and unrealistic in nature. Though there is no empirical
evidence to support this, it is also commonly theorized that the incorporation of ions in the ice
lattice mainly serves to facilitate the formation both ionic and Bjerrum protonic defects through
vacancy or interstitial-like mechanisms.
21,22,23
The ions themselves are assumed to be insoluble
and exit mainly as immobile inclusion or clusters that are securely trapped within the ice lattice.
Therefore, a detailed and realistic atomistic study emphasizing ion dynamics in the ice lattice

4

will significantly address such conjectures regarding the dynamics of ion in ice as well as
providing a further understanding of the true nature of transport processes from the molecular
perspective.
In this investigation, we will also undertake a substantial molecular modeling effort using
the atomistic simulation technique of molecular dynamics (MD) on ice. The MD simulation
systems in our work is comprised of a substantially larger clusters (>6000 molecules) in an effort
to achieve a more realistic interpretation of the simulation data. In order to increase the
likelihood of capturing the unique and critical dynamic processes as they evolve with time, the
simulation time of our study will extend far beyond previous studies in the range of 1 micro-
second ( μ-sec). Detailed analysis of the simulation data will provide substantial insights on the
structural and transport dynamics in ice. We will perform electrical properties calculations on
the simulation experiments in order to offer direct comparison with the laboratory experiments.
Because very few studies currently exist that provides the actual interpretations on the electrical
behavior of ice using atomistic simulations, our combined theoretical and experimental work
represents one of the most comprehensive studies on the physics of ice.
In summary, the goal of this dissertation is to examine the electrical behavior of ice and
doped-ice systems using both EIS and atomistic simulations to develop a comprehensive
understanding of the complex chemical processes involved in the transport properties of ice as it
truly exist in nature. In the following chapters of this investigation, we first describe the
experimental setup and calibration protocols using the EIS measurement technique under
laboratory conditions in Chapter 2. The database of measurement results and characterization
techniques and interpretations through equivalent circuit modeling on ice and ion-doped ice are
discussed in Chapter 3. In Chapter 4, the discussion on the MD simulation development,
experimentation, and results pertaining to transport mechanism in ice is addressed. Finally, a
description of significant scientific findings from our investigation is summarized in the
conclusions of Chapter 5.

6

Chapter 2: Equipment, Experimental, and Calibration

2.1 Background
The fundamental principle of EIS measurements is an application of a non-intrusive,
potentiostatic difference across a material of interest using polarizable contact probes. Since the
relaxation processes are intrinsically unique to different materials, EIS provides a direct and
unambiguous technique to identify and characterize materials by its electrical properties. For
this electrochemical technique to be truly viable for field applications, the contact sense probes
must be constructed of robust material (e.g., electrically conductive and corrosion resistant) and
geometric design in order to eliminate or minimize any measurement artifacts attributable to
electrode polarization and probe degradation. The equations for this electrical measurement
technique for data extraction and analysis are described in Appendix A. By knowing the
geometric factor of the probe configuration, critical frequency dispersive electrical properties
such bulk dielectric constant and bulk conductivity can be readily extracted.

2.2 Experimental
The design and implementation of an electrical system capable of in-situ impedance
measurements on various thermodynamic phases pose formidable challenges. In order to
accommodate a robust electrode configuration suitable for field applications, an electrochemical
properties cups (EPC) was developed to house the electrolyte samples in both liquid and solid
water states. The typical CAD of the EPC is shown in Figure 2.2.1. The standard container
can hold an estimated volume of 30 ml, with holes located at the bottom for insertion of
electrode made from stainless steel screws. The steel screws are also gold-plated for purposes
of mechanical and chemical stability. Since measurements were conducted vertically, o-rings
were used to seal the screw head probes inside the EPC to prevent the liquid sample from
leaking into the electrical hardware. As shown in Figure 2.2.2, actual EPC containers are
constructed from two types of materials: 1) polystyrene for purposes of internal visualizing ice
growth (EPC-1), and 2) Teflon for purposes of long life under extreme low temperatures (EPC-
2 and EPC-3). EPC-1 and EPC-2 designs have identical cell constants, and should produce the
same measurements results. The circuit card beneath the EPC is connected to the 6430B

7

Precision Component Analyzer from Wayne Kerr Electronics, Inc. using 6 feet BNC cables.
Though 4 probes were available, only two terminal connections were used due to expected high
impedance of ice. Multiple probe configuration including parallel plate electrodes from metal
sheets were also adapted to the electrode screws in order to evaluate the effects of electrode
geometry including parallel-plate configuration. For environmental control during
measurements, an environmental chamber (Model ZP-8) from Cincinnati Sub-Zero (CSZ)
Products, Inc. was used to control the temperature. For purposes of experimental consistency
and around-the-clock data collection over a wide range of experimental conditions, a
driver/controller software algorithm was developed to provide automation control of the entire
experimental setup (see Figure 2.2.3). The driver algorithm, written in LabView 7.1, has
capabilities of temperature control of environmental chamber, instrument control of the
impedance system, and experimental control through user-defined soak durations and
experimental description (see Figure 2.2.4).

Figure 2.2.1: CAD of a standard EPC design using Teflon. This is a 2 probe configurations;
similar 4 point configuration were also employed for purposes of adapting other probe
configurations. Unit of length is in inches. Stainless steel screws are used as measurement
probes located the bottom of the EPC container.

8

Figure 2.2.2: Example of an EPC designs with point electrode configuration: polystyrene (EPC-
1) and Teflon (EPC-2 and EPC-3). EPC-1 and EPC-2 designs have the same probe geometry.
EPC-3 is deemed the most robust design suitable for field applications containing recessed
probes.

Figure 2.2.3: Overview picture of the laboratory set-up and equipment list at JPL.

9

Figure 2.2.4: Description LabView program (GUI functionality and source code) for automation
of experimental measurements.

As shown in Figure 2.2.5, electrolyte solutions of KCl, KI, NaCl, MgSO
4
, KOH, and
CaCl
2
were prepared in 50 ml plastic vials. With the exception of MgSO
4
, all salts were used
without additional preparation. MgSO
4
came in hydrated form and was baked for 2 hrs at 80
o
C
prior to use. Typically, an initial concentration of 1 M was prepared for each electrolyte type,
and diluted down to obtain addition concentration levels of 100mM, 10mM, and 1 mM. A 25 ml
sample volume each electrolyte solution was required to perform a single set of EIS
measurements in the EPC container.

10

Figure 2.2.5: Illustration of electrolyte sample preparation.

Measurement conditions were consistent in all experiments. An empty cup calibration
measurement was first performed over the entire frequency range for instrument calibration
purposes and to check for parasitic capacitance. The peak voltage was set to 1 volt amplitude,
with frequency range of 20 Hz to 500 kHz. To obtain high spectra fidelity, 200 data points were
taken per measurement trace; the final value at each frequency was based on the average 3
measurements. Due to current limit considerations, the Wayne Kerr instrument is configured to
measure admittance and phase angle. Typical measurements were taken at 10
o
C decrement
within a temperatures range of 25
o
C to –65
o
C. Each temperature was soaked a minimum of 90
minutes prior to measurement, thus, a single sample requires a total time of approximately 20 hrs
to complete.
Equivalent circuit modeling of experimental results was performed in both ZSimWin
(version 3.10) and Matlab (version 7.8.0.347, R2009a) from Mathworks. ZSimWin, from
Princeton Applied Research, Inc., performs complex non-linear least square (CNLS) fitting of
impedance measurement using Simplex algorithm. ZSimWin is mainly employed in liquid
measurements for evaluation of interfacial behavior by CPE (constant phase element). In order
to optimize accuracy and model flexibility on more complex equivalent circuits, model
algorithms in Matlab were also developed. The more complex equivalent circuits were mainly
employed to analyze the electrical behavior of ice measurements.

11

2.3 Calibration Results & Discussion
Experimental measurements entail both liquid and solid (ice) electrolyte EIS
measurements over a wide thermodynamic range (20 to -65
o
C). The main purpose of liquid
electrolytes is for calibration to ensure electrode configuration provides valid measurement
results at various ionic concentrations. For ice measurements, the ac impedance spectra should
provide frequency dispersion of conductivity and dielectric constant spectra, from which
dielectric relaxation time and activation energies can be ascertained. Analysis and evaluation of
these critical measurement results can potentially provide critical information such as ionic
species and even defect concentration present in ice. Thus, strict calibration protocols must be
developed in order to understand the system limitations and measurement artifacts.

2.3.1 Instrument and Measurement Calibration
All EIS measurements are performed using a benchtop instrument Wayne Kerr 6430B
Precision Component Analyzer capable of extreme current sensitivity down to 100 femto-amps.
We employed extensive calibration protocols to ensure measurement integrity by minimizing
parasitic resistances and capacitances emanating at the electrical leads and connector interfaces.
The first step is calibration at the instrument. As shown in Figure 2.3.1.1, this process involves
the trimming all lead resistance over the entire frequency range due to the four connector lead
associated with the measurement setup. Measurements on the open-circuit and short-circuit
calibration configuration were established using low resistance BNC connectors. Residual lead
resistances and connector capacitances are trimmed-out for all frequencies, and new calibration
coefficients stored in the instrument firmware.

12

Figure 2.3.1.1: Hardware calibration of the measurement setup.

Attempts were also made to account for the remaining parasitic capacitance of the circuit
connector from the EPC probes inside the container were also implemented. This was done to
improve the data quality at the low sensitivity limit especially at high frequencies. This process
mainly entails subtracting the real and imaginary admittance components of the actual sample
measurements from the empty cup measurement at each frequency step. As shown in Figure
2.3.1.2, the empty cup measurement establishes the actual instrument low limit of the
experimental setup. The phase angle is approximately 90
o
over the entire frequency range
indicating good capacitive behavior as expected. The admittance low limit is in the order of 10
-
10
Siemens. Instrument noise is observed at frequencies below 100 Hz, but insignificant over the
critical range of 100 to 500 kHz. Therefore, residual capacitances from the EPC setup appear
minimal, and the sensitivity should be well within the limits of our measurements cases on ices.
Lastly, we compared the measurements results of raw admittance and phase angle data to
that with parasitic capacitance removed from using the empty cup measurements. The
comparison plots along with its corresponding residuals are shown is Figure 2.3.1.3. For the

13

admittance, the residual is less than 6 x 10
-7
Siemens; the residual for the phase angle is less than
1.2
o
. Since pure DI water represents the lower limits of our measurement cases, such small error
in the parasitic capacitance from the EPC setup is deemed insignificant in all our measurement
cases on ices.

Figure 2.3.1.2: Empty Cup measurement of the residual capacitance on calibration.

14

Figure 2.3.1.3: Residuals of calibrated admittance and phase angle on pure de-ionized (DI)
water measurements.

2.3.2 Calibration on Electrode Configuration with Pure DI Water
The measurement results of admittance, phase angle, conductivity, and relative
permittivity spectra using various types of electrode configurations on pure DI water at T = 25
o
C
are shown on Figures 2.3.2.1 to 2.3.2.4, respectively. Comparing the raw admittance spectra,
both plate and point electrodes exhibit similar behavior with a plateau region from 20 Hz to 20
kHz and a sloped region beyond 20 kHz. The long electrodes, however, exhibited some
significant measurement artifacts near the transition region at 20 kHz. This is more prominently
observed in the phase angle spectra. The transition from capacitive region of high phase angle to
the resistive region of low phase angle contains a small secondary plateau, which is often
indicative of film formation or other undesirable processes present.

15

Figure 2.3.2.1: Admittance spectra of Pure DI Water using various types of electrode
configurations at 25
o
C.

Figure 2.3.2.2: Phase Angle spectra of Pure DI Water using various types of electrode
configurations at 25
o
C.

16

The conductivity and relative permittivity of pure DI water from all three electrode
configuration are calculated from corresponding geometric factors also commonly referred to as
the cell constant. These values for the plate, point and long electrodes configuration are
determined to be 8 cm
-1
, 0.8 cm
-1
, 1.5 cm
-1
, respectively. Comparing the conductivity spectra, all
three electrode designs results in a conductivity of approximately 1x10
-6
S/cm for pure DI water.
As expected, the conductivity spectra from the long electrodes exhibited more frequency
dispersion than the measurement results from both plate and point electrodes. The frequency
dispersion is more prominently shown in the relative permittivity spectra. Again, the secondary
plateau is indicated for the long electrode measurement near the 1 kHz to 10 kHz region. This
measurement artifact does not exist for both the plate and point electrode measurements, and
substantiates the use of point electrodes as the preferred electrode configuration in terms of
mechanical robustness and measurement integrity.

Figure 2.3.2.3: Conductivity spectra of Pure DI Water using various types of electrode
configurations at 25
o
C.

17

Figure 2.3.2.4: Relative permittivity spectra of Pure DI Water using various types of electrode
configurations at 25
o
C.

2.3.3 Calibration on Liquid Electrolyte Measurement Results
Further system evaluation was performed using a wide range of liquid KCl electrolyte
concentrations. The admittance and phase angle plot of the entire spectra at each concentration
is shown in Figures 2.3.3.1 and 2.3.3.2, respectively. Both the admittance and phase angle
spectra exhibited no anomalous behavior. As expected, the magnitude of the admittance range
from a high of approximately 10
-1
Siemens to a low of less than 10
-6
Siemens corresponding to 1
M to less than 10
-6
M concentration, respectively. The change in admittance of approximately 1
magnitude matched that of the same change in concentration suggesting both probe functionality
and accuracy of measurement results. The admittance traces also exhibit a capacitive behavior in
the low frequency region. The onset of the capacitive behavior occurs at higher frequencies for
higher concentrations ranging from 100 Hz at 10 milli-Moles (mM) to 10 kHz for 1 M
concentration. The phase angle spectra also indicate resistive or capacitive behavior for 0 to 90
o
,
respectively. At lower frequencies, the capacitive region of the admittance spectra correspond to

18

the relaxation from low to high dispersion in the phase angle. The onset of the relaxation from
the phase angle spectra changes with ionic concentration same as the admittance spectra.
Ostensibly, the low frequency capacitive behavior represents the electrical double-layer
commonly present is all charged surfaces. However, this lower frequency phase angle is
approximately 70
o
, suggesting that the double-layer capacitance does not behave like an ideal
capacitor.

Figure 2.3.3.1: Admittance spectra of KCl electrolytes at 25
o
C. Ionic concentration range is
from pure DI water to 1 M.

19

Figure 2.3.3.2: Phase Angle spectra of KCl electrolytes at 25
o
C. Ionic concentration range is
from pure DI water to 1 M.

The corresponding conductivity and relative permittivity spectra for the same electrolyte
solutions are shown in Figures 2.3.1.3 and 2.3.1.4, respectively. Conductivity spectra are nearly
identical to that of the admittance with differences only in the magnitude due to the geometrical
factor specific to the probe configuration. The conductivity ranges from 10
-1
S/cm at 1 M KCl
concentration to approximately 10
-6
S/cm for pure DI water. Obviously, the conductivity of pure
DI water at 10
-6
S/cm is attributable to residual purities. Further details on the double-layer
capacitance comes from the relatively permittivity spectra. The double-layer capacitance is
comprise of both the water monolayer next to the electrode (known as the stern layer) as well as
a diffuse layer from the concentration gradient emanating from the bulk solution. Also, distance
of closest approach to the electrode by the monolayer of charge species is determined by the size
of the hydrated ions (Barsoukov and ross Macdonald, 2001, pp. 64-66). Therefore, the interface
behaves as two capacitors in series with the diffuse layer so large that it no longer contributes to
the double-layer capacitance and one sees only the water monolayer. It is typically understood
that the double-layer next to the interface exhibits a net dipolar orientation causing the observed

20

dielectric constant well beyond that of bulk solution. From our measurement results on relative
permittivity, this double-layer capacitance for the EPC probe configuration is observed at the low
frequency limit and more prominently in electrolytes with higher ionic concentrations.
Measurement of 1 M ionic concentration provides a relatively constant relative permittivity of
approximately 3 x10
7
starting at a frequency of 10 kHz. Less than 1 M ionic concentration, the
relative permittivity spectra predominantly exhibit a relaxation from a constant value of 80
representing bulk water to that of the double-layer. From this result, the double layer
capacitance can be simply estimated using the equation
l
A
C
o dl
   2.3.3.1
where C
dl
is the double layer capacitance, ε
o
is the permittivity of free space, ε is the measured
relative permittivity, and A/l is the geometrical factor from the ratio of known surface area/length
of the probe configuration. The surface area is determined to be 0.0628 cm
2
and the length
between the probes is 1 cm. Thus, the double-layer capacitance is in the order of 10
-6
F/cm
2
,
which consistent with values reported in literature (Bard and Faulkner, pg. 9) and correspond to a
double-layer thickness in the order of approximately 300 Å for liquids.

21

Figure 2.3.3.3: Conductivity spectra of KCl electrolytes at 25
o
C. Ionic concentration range is
from pure DI water to 1 M.

Figure 2.3.3.4: Relative permittivity spectra of KCl electrolytes at 25
o
C. Ionic concentration
range is from pure DI water to 1 M.

22

Measured results using the calibrated measurement system were also directly compared
with theoretical calculations using known concentrations of KCl liquid electrolytes at T = 25
o
C.
The common equation for calculating the theoretical conductivity (σ) in dilute solutions is
expressed as
i i
i
i
C z F  

 2.3.3.2
where F is the Faraday’s constant of 96,485.34 C/mole, z
i
is the charge number, μ
i
is the
mobility, and C
i
is the concentration of species i. Using equivalent ionic conductivity (λ) of ionic
species i which is defined as
i i
F    2.3.3.3
and readily available from literature (Bard and Faulkner, 2001, pg. 67) is simply calculated at
various ionic concentrations. The summary plot comparing theoretical to measured results at
various ionic KCl concentrations is shown in Figure 2.3.3.5. The measured results compare
reasonably well to theoretical values especially at ionic concentrations greater than 10
-1
mM
which indicate high solubility with chemical activity of near unity. At ionic concentrations
lower than 10
-5
mM, higher deviations from theoretical were observed. This is a consequence of
a combination of both residual impurities and the natural ionic concentration of approximately
10
-4
mM in the form of H
+
and OH
-
ions present in the DI water (Fletcher, 1970, pg. 149).

23

Figure 2.3.3.5: Comparison of theoretical vs. measured conductivity for various concentrations
of KCl electrolytes at T = 25
o
C. The equivalent conductivities used here are: λ (K
+
) = 73.25
cm
2
-S/-mole, λ (Cl
-
) = 73.25 cm
2
-S/mole. Measured conductivities were extracted at 100 kHz
frequency.

24

2.4 Conclusions
The design apparatus capable of in-situ EIS measurements on ice using polarizable
electrode probes is highly successful in terms of sensitivity range and accuracy. Calibration of
the system to remove both parasitic resistances and capacitances performed at the lead and the
circuit connection to the baseplate of the EPC. With or without calibration at the connection, the
residual difference of DI water in admittance is less than 10
-7
S while the residual difference in
phase angle is less than 1.2 degrees. Thus, residual capacitance at the circuit connector is
insignificant; compensation of parasitic resistance and capacitance of the connector leads remain
dominate factors in measurement accuracy. Comparison of various electrode configurations, the
most robust design in terms of simplicity, accuracy, and mechanical integrity was the point
electrodes. The long electrodes performed the worst due to the significant frequency dispersion
observed in the bulk measurement region. Comparison results using various concentrations of
KCl electrolytes ranging from 10
-4
mM to 1 M also exhibited excellent calibration accuracy. A
change in one magnitude of ionic concentration corresponds directly to the same change in bulk
admittance and conductivity. From the phase angle spectra, the onset of the electrical double-
layer was detected, and changes with ionic concentration. The relative permittivity of the
double-layer is determined to be approximately 10
7
, which corresponds to double-layer
capacitance of 10
-6
F/cm
2
with a double-layer thickness of approximately few hundred
angstroms. Our results are within expected values reported in literature. Lastly, measured
versus theoretical values of bulk conductivities indicate high solubility with a chemical activity
of near unity, with deviation occurring at the lowest ionic concentration. Such deviation is
expected and is mainly attributable to impurities present. Combined with strict calibration
protocols to ensure optimum measurement sensitivity and integrity, our experimental system
exhibited full functionality and capabilities of in-situ ice measurements overall.

25

Chapter 3: Electrical Spectroscopy Measurements of Ice

3.1 Theory
Though the ice structure has many structural forms, the most common is ice Ih. The ice
Ih structure has four nearest neighbors at the corners of a regular tetrahedron. The hydrogen
atoms, are covalently bonded to the nearest oxygen to from H
2
O molecules, and these molecules
are linked to one another by hydrogen bonds; each molecule offering its hydrogen atoms to two
other molecules and accepting hydrogen bonds from another two. There are two possible
hydrogen sites on each bond with four of these sites adjacent to each oxygen. The disorder of
the hydrogen atoms over these sites satisfies the two ice rules: 1) there are two hydrogen atoms
adjacent to each oxygen, and 2) there is only one hydrogen atom per bond.
The classical theory on ice transport properties is focused on crystalline defect since most
ionic species are virtually insoluble in ice, and when present they mainly exist as inclusions or
clusters. However, the electrical properties of ice are thought be to be very sensitive to even
small concentrations of certain types of ionic species or impurities from studies on HCl, NH
3
,
KOH, and most importantly HF. It is from these specific case studies that formed the basis for
the classical understanding of the electrical properties of doped ice today. This basic
understanding mainly states that ionic species incorporated in the hydrogen bonded network
facilitates the generation point defects which can cause local reorientation changes of the ice
lattice. Such defects must locally violate the ice rules thereby reorienting the molecules from an
ideal structure to another along a defined path. Thus, this classical interpretation emphasizes the
type, concentration, and mobility of such defects as the main contributor to the polarization and
conduction processes, while the ionic species themselves do not exhibit appreciable transport if
solvated in the lattice.
To distinguish them from conventional vacancies and interstitials defects, the point
defects thought to be responsible for the electrical properties in the ice lattice are referred to as
potonic point defects. There are two pairs of protonic point defects that can be formed in a
perfect ice structure, making four types of defects in total. The first pair is known as the Bjerrum
defect pair. As shown in Figure 3.1.1, these two types of protonic point defects are the results of

26

a re-orientational turn of a single water molecule producing on bond with no protons and another
with two protons towards one another. The empty bond is an L-defect, and the bond with two
protons is the D-defect. Such ‘orientational faults’ are considered unstable and typically
inducing additional molecular turns and further separation of the defective bonds as dictated by
local lattice energetics. The second pair of protonic point defects is the H
3
O
+
and OH
-
ionic
defects. This defect pair is formed by the transfer of proton from one molecule to the neighboring
molecules, and separated by successive jumps of protons from one end of the hydrogen bond to
the other. The common ionization reaction
 
  OH O H O H
3 2
2
is a typical process in liquid water, but in ice the ions do not move as complete units. The
oxygen cannot migrate from one site to another, and the motion of the proton along a bond
transfer the state of ionization from one molecule to another. Therefore, the polarization and the
electrical conduction mechanisms observed in electrical measurements are the direct
consequence of the interactions of both Bjerrum and ionic protonic defects in the ice lattice.

Figure 3.1.1: Illustration of L,D protonic Bjerrum in the ice lattice. The re-orientation turn of a
single water molecule forms both L,D-defect.

27

3.1.1 The Debye Theory for Relaxation in Ice
The electrical properties of materials as determined by the alternating current admittance
(reciprocal of impedance) exhibiting a single relaxation is commonly referred to as the Debye
relaxation (Petrenko and Whitworth 1999, pp. 64-68). The Debye relaxation mechanism in ice is
defined by a slow process that can be accurately described by the general polarization equation
for an electrical response
  P E
dt
dP
s
D
   

0
1
3.1.1.1
where P is the polarization, χ
s
is the susceptibility, E is the applied electric field, ε
o
is the
permittivity of free space equal to 8.85x10
-14
F/cm, and τ
D
is the Debye relaxation time constant.
From equation 3.1.1.1, the polarization is dependent of two parameters: the strength of
susceptibility, χ
s
, and its time constant τ
D
. The susceptibility, χ
s
, in ice defines the polarization
strength and can also be described by the equation

    
s s
3.1.1.2
where ε
s
is the low frequency relative permittivity (also known as the “static” value) and ε
∞
is the
high frequency relative permittivity (also known as the “infinite” value).
For a sinusoidal alternating filed of angular frequency (ω), the electric field is now
represented by a complex quantity
t j
o
e E E

 3.1.1.3
where E
o
is the amplitude, ω = 2πf, and 1   j . The polarization response, therefore, is
t j
o
e P P

 3.1.1.4
where P is also a complex value. Simple mathematics leads to a complex frequency-dependent
susceptibility

28

 
D
s
o o
j E
P



 

 
1
. 3.1.1.5
Separating the χ(ω) into the real and imaginary components
            " ' j   , 3.1.1.6
we get the following
 
2 2
1
'
D
s
 

 

 3.1.1.7
 
2 2
1
'
D
s D
 
 
 

 . 3.1.1.8
Ice can also be represented as a resistor. In this case, we employ current density (J)
t j
o
e P j
dt
dP
J

   . 3.1.1.9
The complex conductivity is then
        
o
j
E
J
  . 3.1.1.10
The real and imaginary parts of the complex conductivity is
            " ' j   , 3.1.1.11
where
 
2 2
2
1
'
D
s o D
 
   
 

 3.1.1.12
 
2 2
1
"
D
s o
 
 
 

 . 3.1.1.13
A more general derivation of the frequency-dependent Debye equations from first principles is in
Appendix B.

29

3.1.2 Polarization and Resonance in Ice from a Molecular Perspective
At the most basic atomic or molecular perspective, the application of an electric field
does not induce instantaneous polarization because charge particles have inertia. In the case of
atomic and ionic polarization, the behavior of electrons and ions, to a first approximation, are
analogous to linear springs bound to equilibrium positions. Thus, the restoring force is
proportional to the displacement x from their equilibrium position. Since the mechanisms for
polarization is accompanied by energy dissipation, a damping factor is typically included in the
equation of motion. If the damping (resistive) force is assumed proportional to the velocity of
the moving charged particle of mass, m, displaced by distance, z, from its equilibrium, the full
equation of motion can be described as

m
E q
z
dx
dz
dt
z d
n

  
2
2
2
2  
3.1.2.1
where 2α is a damping factor, ω
n
is the system natural resonant angular frequency, and E’ is the
local electric field as given by the Mosotti approximation:

0
3 
P
E E    3.1.2.2
where P is the polarization charge density in C/cm
2
. The external or driving steady-state
sinusoidal electric field, E, is given by the same equation 3.1.2.3.
The polarization, P, is described by
eNz P  3.1.2.3
where e = 1.6x10
-19
C and N is the number of oscillators per unit volume.
Substituting the above relations into the harmonic oscillator differential equation yields

2
2
0 2
2
d P dP
P kE
dt dx
   
3.1.2.4
where the k-coefficient is:
m N e k /
2
 3.1.2.5

30

and the resonant angular frequency is:
m
N q
n
0
2
2 2
0
3 
   
3.1.2.6
indicating that the resonant frequency is less than the natural frequency.
The solution to the polarization differential equation is:
22
0
2
kE
P
j    


3.1.2.7
which implies there is a phase difference between the P and E vectors.
Next, we use the polarization equation to model the orientation polarization found in
water/ice. At high frequencies, the polarization is due to the sum of the electronic and atomic
polarizations and is given by P
∞
= ε
0
ε
∞
E. The water/ice orientation polarization adds to the
electronic and atomic polarizations and is given by:
 E P

    
0
3.1.2.8
Since ε = ε
s
at ω = 0, combining the previous two equations at ω = 0, leads to the intermediate
result
 

    

s o
o
k
2
3.1.2.9
Combining the above three equations and the Debye time constant, τ
D
= 2α/ω
0
2
, leads to an
expression for the relative permittivity

 
D o
s
j  
 
 
 

 

 2 2
/ 1
* 3.1.2.10
Separating the real and imaginary components of the complex relative permittivity
" ' *    j   , 3.1.2.11
you get

31

   
 
2 2
2
2 2
2 2
/ 1
/ 1
'
D o
o s
   
   
 
 
 
 


3.1.2.12
 
 
2 2
2
2 2
/ 1
"
D o
s D
   
  

 



. 3.1.2.13

3.1.3 Conduction in ice
Conduction of species i (σ
i
) is constant at constant temperature and is defined by
i i i i
Q N    3.1.3.1
where N
i
is its concentration, Q
i
is the charge, and μ
i
is the mobility of species i. Thus, the
fundamental theory for ice conduction is similar to all crystalline materials and depends on the
presences of vacant sites into which defects and/or ions can move. The “extended” form of the
Nernst-Einstein equation as described by
T k
Q N
D
b
i i
i
i
2


3.1.3.2
links the addition of an electric field to the transport of species i using D
i
as the diffusion
coefficient, T is the temperature, and k
b
is the Boltzmann constant of 8.62x10
-5
eV/K. The
mobility (μ
i
) of a defect i is defined by its drift velocity in the electric field, and thus, are
normally expected to vary with temperature according to an Arrhenius relation of the form








 
T k
E
T
b
im
i
exp
1
 3.1.3.3
where E
im
is the free energy of activation for motion. Additional factors that typically contribute
to mobility are charge and size of the conducting species and lattice structure. A highly charged
species will polarized the ions of opposite charge as it migrates pass by, inducing an increase in
the energy barrier that further inhibits the change of site. Larger species tends to also hinder in a
similar way by the interaction of its outer electrons with nearby species it must pass in order to
reach a new site (Moulson and Herbert, 1992, pp. 43-46). Some lattice structures may provide
space for movement through channels to facilitate motion of certain charged species.

32

The conductivity of the crystalline material can be described in terms of defect
concentration in the crystal structure. From Ohm’s law, conductivity is simply the
proportionality constant between the current density (j) and the applied electric field (E) such as
E j   3.1.3.4
From Equation 3.1.3.3, the expression for current density can alternatively be expressed by








 
T k
E
E
T
A
N j
b
im
LD
exp 3.1.3.5
where N
LD
is the Bjerrum defect concentration and A is constant describing the vibrational state
of the crystal. The temperature dependent Bjerrum defect equilibrium concentration (N
LD
) in the
lattice is often described by the form








 
T k
E
N N
b
LD
LD
exp 3.1.3.6
where N is the number of molecules per unit volume available in the lattice and E
LD
is the
associated activation energy accompanying the defect formation. Thus, conductivity equation
can also be re-written in an Arrhenius form








 
T k
E
b
t
o
exp   3.1.3.7
where σ
o
is a temperature independent pre-exponential coefficient and E
t
is the total energy
associated with defect formation and migration induced from the electric field as described by
equations 3.1.3.3 and 3.1.3.6.

33

3.2 Literature Review on Ice Electrical Properties
The electrical properties of ice have a long history of distinguished investigators
beginning with Debye (Debye, 1929). With simple models of spherical dipoles, Debye was able
to deduce the mechanism of rotating or re-orienting H
2
O molecules under the combined
influence of an applied electrical field and the inner viscosity of ice to account for the observed
relaxation processes as a function of frequency. Debye also rationalized that only a few
molecules are free to reorient (~1 in 5 million) to retain stability of the ice lattice, but was not
able to provide an actual value for the number of rotating molecules in ice. The Debye Theory
also did not provide an accurate description of the anomalously high mobility observed for H
+

and OH
-
ions in water or ice. The study of Bernal and Fowler offered the first molecular model
for potential energy surface (PES) in the development of the quantum mechanical theory as a
mechanism for ionic mobility.
24
This theory attributes the unusually high mobility of H
+
and
OH
-
ions observed in water to a combination of orientation of dipole molecules along with
proton hopping from molecule to molecule. Transport of electrical current in ice was not
addressed with significant detail until the work of Bjerrum in 1952.
25
In this seminal study,
Bjerrum employed both theories of Debye and Bernal and Fowler to describe the ability of ice to
conduct current through the existence of thermally induced movements of electric dipoles
initiated at the defects of the ice lattice. The two types of defects are: (1) orientational, where
two water molecules are incorrectly oriented resulting in either two (D-defect) or no (L-defect)
protons between adjacent oxygen atoms, and (2) ionic, where H
3
O
+
and OH
-
ions exist in the
lattice instead of H
2
O molecules. The passage of current is the result of migrating H
3
O
+
and OH
-

ions in an electric field counteracted by movement of dipolar water molecules in the opposite
direction until charge recombination site in the lattice is reached. The work of Bjerrum offered
the first complete mechanism for molecular ice dynamics in support of the phenomenological
electrical properties, which laid the foundation for modern interpretations of the importance of
lattice defects in pure and doped ice.
Dielectric properties correspond directly to the extent of polarization due to an applied
electric field. Molecular polarization in ice is related to the degree of proton disorder within its
lattice. The dielectric properties of ice were first explored by Pauling in 1935.
26
Pauling
hypothesized that water molecules in ice form a hydrogen bond network with four neighbors

34

with each hydrogen bonded to an oxygen and a hydrogen bonded to another. The many different
arrangements of the hydrogen were substantiated by the residual entropy of ice at low
temperatures observed experimentally. Kirkwood was the first to provide modern interpretation
between the measured dielectric constant and its molecular constituents in 1936.
27
Based on the
work of Lorentz who produced the first estimate of the average local field in a molecule,
Kirkwood established the statistical calculation of dielectric polarization from a molecular point
of view in terms of the average induced moments from uniform fields acting in its interior. This
relationship became to be known as the Kirkwood-Fröhlich equation for infinite systems of its
own dielectric constant such as water and ice (Fröhlich, 1949).

More simplified treatment for
computing static dielectric constant of ice employed a spherical region inside the specimen large
enough for the material outside to be treated macroscopically.
28
The integrated result based on
the configuration energy of the system is a statistical average of electric moment of a sphere or
group of spheres. Neumann, later, provided a full theoretical description

for dielectric behavior
with the addition of an applied electric field, thereby providing formulations on computing
dielectric constant due to the combined polarization of both external and local fields.
29

Early experimental studies mainly focused on the elastic and dielectric behaviors
observed from DC as well as AC electrical responses. However, reliable data was often difficult
to achieve due to impurities as well as known ambiguity or errors associated with anisotropy and
crystalline orientation in reference to the applied electric field.
30, 31
To mitigate such measurement
artifacts, almost all early investigations employ a parallel plate electrode configuration. Autry
and Cole provided the first accurate study on the dielectric behavior of an ac response and
determined that the relative permittivity of ice at high and low frequencies were approximately 4
and 91.5 at 0
o
C.
32
The static permittivity was found to increase with decreasing temperature as
expected due to changes in degree of proton disorder. Later, Granicher was credited with
deducing the similarities between energies and jump times that define diffusion and dielectric
relaxation process.
33
This pivotal study enabled Haas to correlate activation energies to
diffusive mechanism due to lattice imperfections in ice and approximated a defect concentration
of 10
13
cm
-3
at 0
o
C.
34
Onsager and Dupuis were the first to interpret the electrical properties of
ice on the basis of electrical moments from interactions between complexes of adjacent defects
(Onsager and Dupuis, 1962, pg. 27). Onsager and Runnels further examined the underlying
mechanisms for diffusion and dielectric relaxation by cross comparison of several rate processes

35

due to ice imperfections, and concluded that interstitial type migration of an ‘intact’ water
molecule is the most likely mechanism for diffusion of defects in ice lattice.
35
Though never
experimentally conformed, this study specifically outlined the likely composite configurations
resulting from L,D defects originating from interstitial molecules. Recently, a mechanism for
proton transport originating from interstitial diffusion of intact water molecules was
experimentally supported by work with NMR.
36
In relating the dielectric constant of ice to
defects concentration, Nagle the first to surmise that even a small number of ice defects play a
crucial role in the dielectric constant of ice. He also suggested that ice defects only affect the
rate at which ice responds to an applied field or relaxation phenomena but not the final
polarization.
37
The result of this work offered a quantification of effective partial charges of
Bjerrum and ionic defects in ice of 0.38e and 0.62e, respectively, and is still generally accepted
today. The critical work of von Hippel strived for further clarification of ice polarization and
conduction by means of molecular interpretations of single crystalline ice.
38
Based on Jaccard
theory which takes into account both L, D and ionic defects, the observed dispersion in
polarization and conductivity is determined by “majority carriers” (defect for which
“concentration x mobility” is largest) which could be analyzed through relaxation times.
39
At
high frequencies, the majority carriers are Bjerrum defects in motion due to thermal activation;
ionic defects are dominant at low frequencies due to its intrinsically higher mobility. However,
the effort to offer molecular insights to measurement results turned out to be greatly complicated
by experimental conditions to grow and to measure ice which could account for the effects of
anisotropy and, possibly, eliminate the distribution of relaxation times arising from Maxwell-
Wagner interfacial polarization.
40, 41

Studies on doped ice also provided further insights into the complexities involved in the
mechanisms of defect generation and transport under more realistic conditions and lower
temperatures. Gough and Davidson were the first to attribute the observed sensitivity of
electrical properties to impurities resulting from a buildup of extrinsic Bjerrum defects observed
at low temperatures.
42
This increase in Bjerrum defects greatly outnumbers the intrinsic defects
associated with the ideal ice lattice and is indicated by a dramatic decrease in relaxation times.
Much of today’s understanding of the electrical properties of ice originated from the seminal
study by Camplin with HF-doped ice was the first comprehensive attempt to apply the theoretical
models of current conduction as it relates to defect concentration.
43
With known impurities

36

concentration of HF, the concentration of L defects increased due to substitution of H
2
O
molecules in ice, thereby, creating proton free bonds in the lattice. Thus, the observed low and
high frequency conductivity can be interpreted through variances in concentration, charge,
mobility, and activation energy of defect formation associated with Bjerrum and ionic defects.
This study further validated previous works of Jaccard and von Hippel in characterization of
electrical properties through the interpretation of molecular defect behavior. Working with HCl
and KOH-doped ice, Kawada et al. more recently added further details on the mechanisms of
dielectric relaxation from activation energies of shortened relaxation times involving the L and
ionic defects associated Bjerrum defect complexes with lattice imperfections such as vacancies
and interstitials at wider temperature ranges (Kawada et al, 1989 and 1997). In addition to
further supporting the critical role of Bjerrum defects cited in earlier studies, it was also found
that partial phase changes were likely to occur below 110
o
K due to constraining of the molecular
rotational mode, resulting in possible existence of inhomogeneous phases at 75
o
K. Johari and
Whalley confirmed the inference of the orientational ordering by the measurement of dielectric
relaxation spectra in the 0.05 Hz-300 kHz range, and extracted the relaxation time (τ) of 1.5 sec
with a corresponding activation energy of 23.4 kJ/mole.
44
Thus, this observation provided some
empirical evidence in support of the type of defects that are mobile and consequently play a
dominate role in the electrical properties in ice.

37

3.3 In-Situ Ice Observations
Prior to actual EIS measurements on pure and more importantly doped ice, we performed
a thorough visual inspection on the ice crystal structure to gain some initial insights to correlate
with measurement behavior. The ice is slowly grown in-situ at 5
o
C decrements down -20
o
C
using the clear polystyrene EPC design. After freezing, we took close-up pictures of the ice at
several angles inside the EPC to deduce the in-situ freezing mechanism. We also took pictures
of the ice over a light source to enhance visibility of the various phases between single crystal,
ice fractures, and precipitated salt through differences in reflectance. Figure 3.3.1 illustrates the
baseline visual characteristics of the freezing process of pure ice grown from DI water. From
such pictures, it is understood that ice nucleation uniformly occurs at the EPC outer wall and
propagates toward the center causing ice figuring and/or fractures as a consequence of material
insolubility in the ice lattice. The solid cluster observed in the middle is attributable to residual
impurities in the DI water. In addition, it is observed that the ice at the bottom where the probes
are located is relatively free of fractures which suggest the desirable pure single crystal ice
measurements as oppose to poly-crystalline type near the center of the EPC.
From a basic perspective, ion distribution within the ice lattice should follow the
nucleation pathway thereby significantly impacting the defect concentration and distribution.
This is definitively observed in ice doped with various concentrations of K
+
and Cl
-
ions.
Compared with pure ice as shown in Figure 3.3.2, ice doped with higher ionic concentrations
exhibited less opaqueness towards the EPC center attributable to higher salt precipitation and/or
ejection from the ice lattice during nucleation. The salt cluster from 100 μM KCl solution
(Figure 2.3.3) has radius 2-3 times less that than of that from 100 M KCl solution (Figure 3.3.5),
and nearly covers the entire radius of the EPC. With 1 M KCl solution as shown in Figure 3.3.6,
the salt cluster appears uniformly distributed in the ice most likely due to ionic saturation. Thus,
higher ionic concentrations may results in amorphous-like ice similar to liquid water.
Upon placing light source at the bottom of the EPC, the relative transparency to visible
light and coloration through the ice varies substantially with ionic concentrations. On pure ice,
the light is relatively bright (Figure 3.3.2b). For 100 μM and 10 mM concentration, the color
exhibited a darker color towards light red (Figure 3.3.3b and 3.3.4b). The lowest transparency of
light is exhibited by 100 mM and 1 M KCl concentration (Figure 3.3.5b and 3.3.6b) with the
color of dark blue. Solid color is exhibited in all ice sample near the center of the EPC, and is

38

attributable to salt cluster ejected from the lattice during nucleation. The mechanism for such
observed optical differences is not addressed in this study, but from classical theories it is
attributable to ionic absorption and ejection processes in the ice lattice which can potentially
have a substantial impact on phenomenological interpretation from EIS measurements.

Figure 3.3.1: Ice grown from pure DI water (side view).

Probes footprint

39

Figure 3.3.2: (Left) Ice from pure DI water (top view) without light source underneath. (Right)
Ice from pure DI water (top view) with light source underneath.

Figure 3.3.3: (Left) Ice doped with 100 μM KCl (top view) without light source underneath.
(Right) Ice doped with 100 μM KCl (top view) with light source underneath.

40

Figure 3.3.4: (Left) Ice doped with 10 mM KCl (top view) without light source underneath.
(Right) Ice doped with 10 mM KCl (top view) with light source underneath.

Figure 3.3.5: (Left) Ice doped with 100 mM KCl (top view) without light source underneath.
(Right) Ice doped with 100 mM KCl (top view) with light source underneath.

41

Figure 3.3.6: (Left) Ice doped with 1 M KCl (top view) without light source underneath. (Right)
Ice doped with 1 M KCl (top view) with light source underneath.

42

3.4 Measurement Results on Ice
EIS measurements of doped ice were performed using the experimental setup outlined in
Chapter 2 from liquid electrolytes containing KCl, CaCl
2
, KI, KOH, NaCl, and MgSO
4
and ionic
concentration ranging from 1 μM to 1 M. Measurements of pure ice from DI water using both
plate and point electrode configurations were also performed to establish baseline in terms of
measurement accuracy. From the raw admittance and phase angle measurements, calculations of
loss tangent, relative permittivity, and conductivity were calculated from the known geometric
factor as determined from calibration protocols. As a general practice, Nyquist-type plots on the
real and the imaginary components of admittance were employed to establish the measurement
accuracy and integrity of each ice sample. In addition, loss tangent, conductivity, relative
permittivity plots over a wide thermodynamic range were also summarized for comparison
purposes.
Plot summaries for pure ice using both parallel plate and point electrodes are shown in
Figures 3.4.1a-e and 3.4.2a-e, respectively. The magnitude and frequency dispersion of the
admittance spectra of pure ice using parallel and point electrodes are almost identical. The
frequency of the Debye relaxation for both electrode configurations decreases with decreasing
and ranging from 10 kHz at -15
o
C to 100 Hz down to -65
o
C. The point electrode design does
exhibit signatures of a second relaxation at high frequencies less than 100 Hz. This is an
indication of a liquid water layer from thermal heating due to the applied potential difference
(p.d.). From the Nyquist plots, the point electrodes also exhibited an inductance loop. The
observed inductance is only present in ice, not liquid water. This unique measurement signature
is due to the ice resonance at onset of the Debye relaxation, and is attributable to the radial
orientation of the applied electric field. This magnitude of the resonance (called Debye
resonance going forward) is also thermodynamic. Unlike the Debye relaxation, it increases with
decreasing temperature. The Debye resonance is also exhibited in the relative permittivity
spectra as indicated by the “dip” at the onset of the Debye relaxation from the equilibrium ε
∞
= 3
at high frequencies to ε
s
=80 at low frequencies. The conductivity of ice for both electrode
configurations ranges from approximately 10
-7
to 10
-9
S/cm from -15
o
C to -65
o
C, respectively.
In the high frequency range, the conductivity is substantially more temperature dependent than in
the low temperature range less than 1 kHz. The classical interpretation for this behavior is that

43

L,D-Bjerrum protonic defects is thermodynamic but not ionic protonic defects. This is in
agreement with our measurement results as Bjerrum defects should dominate in the high
frequencies while ionic defects tend to dominate at low the frequencies. In the region dominated
by Bjerrum defect, however, results from the point electrodes exhibited a small negative slope in
the high frequency range. The loss tangent peak is an indication of the Debye relaxation, and it
increases with decrease in temperature from 4 at 15 for the point electrodes while the peaks of
the loss tangent spectra from the parallel plate configurations are approximately 3 and
temperature independent. The added signatures from the point electrodes are additional
detection parameters of ice should provide further insight on the electrical properties of ice and
doped ice as described in subsequent sections.

3.4.1 Ice Doped with KCl
Ice grown from KCl electrolytes represents the bulk of our doped ice measurements as
such ions are known to incorporate readily into the ice lattice leading to significant changes to
the ice relaxation processes. The electrolyte concentrations range from 1 μM to 1 M. At the
lowest concentration of 1 μM KCl is shown in Figures 3.4.3a-e. The Nyquist plots over the
entire temperature range are similar to that of pure ice including the inductance loop at low
frequencies. As expected, the admittance spectra show a small increase over pure ice at higher
temperatures in low frequencies attributing to mass-transfer effects from an ionic concentration
gradient. Similar to pure ice, the conductivity spectra decrease dramatically with decreasing
temperature in the high frequency range but not as much in the low frequency range. Both the
relative permittivity and loss tangent spectra are similar to pure ice as well. The equilibrium
polarization states at low (ε
s
) and high (ε
∞
) frequency ranges are also similar at approximately 80
and 3, respectively. The loss tangent peaks exhibited the same magnitude ranging from 4 to 10.
From the peak loss tangent, the Debye relaxation ranges thermodynamically from 10 kHz to 100
Hz like pure ice. Therefore, ice doped with small concentration of KCl as does not appear to
appreciably change the ice lattice structure as interpreted from EIS measurements.
Dramatic differences in measurement results in ice were observed with dopant KCl
concentrations of 100 μM and 10 mM as shown by Figures 3.4.4a-e to 3.4.6a-e. As expected, the
Debye relaxation frequency increased with ionic concentration. However, the ice doped with

44

100 μM KCl exhibited the most liquid-like behavior over the entire temperature range down to -
65
o
C. Though Nyquist plot contains an inductance loop like ice, the low frequency ε
s
value is
substantially higher than that of ice at approximately 1000. As the dopant concentration
increased to 10 mM, the ε
s
decreased to 300. The Debye relaxation as indicated by the peak loss
tangent also varied substantially from approximately 20 to 40 for dopant concentration of 1 mM
and 10 mM, respectively. This suggests that the distribution of K+ and Cl
-
ions and their
molecular interactions during the freezing process play a critical role in ion absorption into the
ice lattice at the low frequencies. Such behavior indicates that ions are accumulating at the
electrodes rather than trapped in ice cavities, thus inhibiting ice freezing at the electrodes. From
our measurement results, optimal KCl dopant concentration is 100 μM where the highest ion
absorption in the ice lattice is observed resulting in the most water-like behavior.
At the highest KCl dopant concentrations, measurement results indicated ice conditions
similar to pure ice crystalline structure. As shown in Figures 3.4.7a-e and 3.4.8a-e, spectra of
admittance, relative permittivity, conductivity, and loss tangent exhibited nearly the same
behavior as that observed in pure ice. The ε
s
value agrees with what is expected at
approximately 90. The ice containing 1 M KCl dopant concentration exhibits ε
s
beyond 100,
which can be attributed to excess charge near the electrode interface as the ions are ejected from
the ice lattice during freezing. The peak loss tangent is relatively constant with values in the
range of 4 to 7 and decreases in frequency with decrease in temperature. The Debye relaxation
frequency are the same for both dopant concentrations of 100 mM and 1 M and increased
dramatically for 100 kHz at -25
o
C to 1 kHz at -65
o
C, and is attributable to increases in protonic
defect concentration. Another significant difference is the high frequency conductivity which is
approximately 1 order of magnitude higher than that of pure ice at 10
-6
S/cm. Thus, ion
absorption tends to achieve saturation at high dopant concentrations with the majority of K
+
and
Cl
-
ions get ejected from the lattice during freezing and recombine into the observed salt clusters.
Though the classical interpretation on the increase in conductivity is theoretically attributable to
the increase of Bjerrum defects
45
, it is unobtainable by EIS as to whether such high
concentrations of K
+
and Cl
-
ions are truly immobile within the ice lattice and its relative
contributions to the ice transport as measured by conductivity.

45

3.4.2 Ice Doped with KOH
The lowest dopant concentration of KOH of 1 μM also exhibited ice characteristics as
shown in Figures 3.4.9a-e. The Nyquist plots have the inductive loop indicating the Debye
resonance as also detected by the dip in the relative permittivity spectra. The ε
s
and ε
∞
values at
the equilibrium polarization are approximately 100 and 3, respectively. The maximum high
frequency conductivity at -15
o
C is also indicative of ice at approximately 10
-7
S/cm, and
decreases down to 10
-9
S/cm at -65
o
C. The peak loss tangent is approximately 7 and is
relatively independent with temperature. However, the Debye relaxation frequency is similar to
pure ice ranging from 10 kHz at -15
o
C to 100 Hz at -65
o
C. The 10 mM KOH dopant
concentrations still exhibited ice characteristics as well as shown in Figures 3.4.10a-e. The
inductance loop is present in the Nyquist plots, and the ε
s
is 100 along with 3 for ε
∞
. Though the
conductivity is in the range of 10
-7
S/cm at the high frequencies, the low frequency conductivity
is nearly 1 order of magnitude higher at 10
-8
S/cm which is indicative of ion absorption ice.
Water-like behavior is exhibited at higher dopant concentrations as a consequence of excess K
+

and OH
-
ions prohibiting the freezing process. From the 100 mM KOH measurements as shown
in Figures 3.4.11a-e, Debye relaxation and resonance is also observed. Though the high
frequency ε
∞
value of 3 indicates the presence of ice, the relative permittivity at high frequencies
is in the order of 1000 and beyond at temperatures down to -65
o
C. The measured conductivity
values is substantially higher than that of ice as high 10
-4
S/cm at -25
o
C. Thus, we can conclude
that the solubility of KOH in ice substantially higher causing with liquid-like behavior at
temperature down to -65
o
C.
3.4.3 Ice Doped with CaCl
2

The majority of measurements on ice doped with CaCl
2
also exhibited liquid-like
behavior. Measurements on 1 mM dopant concentration exhibited ice characteristics in the high
frequency range (see Figures 3.4.12a-e). The Debye relaxation is detected by the large peak loss
tangent greater than 20 at frequencies in the order of 10 kHz. It is also observed that the Debye
frequency increased similarly to the ice doped with KCl at 100 kHz near the concentration of 100
mM and beyond. The high frequency ε
∞
value of 3 is also a characteristic of ice, but two

46

relaxation processes at lower frequencies are observed. The two relaxations are expected to be
the nominal ice Debye relaxation and the Maxwell polarization next to the electrodes (also
known as space charge effects). In the Nyquist plots, the Maxwell polarization can be indicated
by the elliptically shape due to the merging of two relaxation semi-circles in the low frequencies.
The most prominent indications of Maxwell polarization which can be distinguishable from
Debye relaxation are the 10 mM dopant measurements as shown in Figures 3.4.13a-e.
Measurements containing Maxwell polarization typically do not reach equilibrium down to the
frequency limit of 20 Hz, with values substantially beyond the characteristic ice ε
∞
value of
approximately 100. The conductivity of CaCl
2
doped ice is almost 1 order of magnitude higher
at 10
-6
S/cm At dopant concentrations of 100 mM and beyond, the conductivity is again 3-4
orders of magnitude higher than that of pure ice (see Figures 3.4.14a-e) which suggests liquid-
like behavior with substantial contributions from free ion motion. It is observed that freezing
does not occur in electrolytes with greater than 10 mM CaCl
2
down to temperatures of -55
o
C
attributable to freezing point depression from high solubility of Ca
+2
and Cl
-
ions in ice.
3.4.4 Ice Doped with KI, NaCl, and MgSO
4

To further investigate the dependency on ionic size and solubility near such saturated
conditions, measurements on 100 mM dopant concentrations of KI, NaCl, and MgSO
4
were also
conducted over the same thermodynamic range. Comparative plots are shown in Figures
4.4.15a-e to 4.4.17a-e. As expected, ice was observed at different sub-zero temperatures
depending on the ionic type. Measurements of both KI and NaCl dopant types exhibited ice
signatures at -35
o
C, and for MgSO
4
at -25
o
C. All the ice signatures are very distinct with no
significant Maxwell polarization present. All Nyquist plots exhibited the thermodynamically
inductance loop which corresponds to the Debye resonance. The ε
s
and ε
∞
values for all dopant
types are approximately 100 and 3, respectively. In frequency range dominated by Bjerrum
protonic defects, the conductivity of MgSO
4
is more than 1 order of magnitude lower at 10
-7

S/cm compared to KI and NaCl doped ice at the same temperature of -35
o
C. The low frequency
range dominated by ionic protonic defects for MgSO
4
is also nearly identical to that of pure ice,
while over the same frequency range KI and NaCl doped ice converged throughout all
temperatures. From the loss tangent, Debye relaxation frequency for the ice doped with MgSO
4

is very similar to pure ice, but the magnitude of the loss tangent is significantly higher at

47

approximately 15. The loss tangent peaks for both KI and NaCl doped ice are less 10. In
addition, the Debye relaxation frequencies for KI and NaCl doped ice have increased similar to
KCl doped ice at approximately 100 kHz. Overall, these signatures strongly suggests a poor
solubility of MgSO
4
salt in the ice lattice as compared to salts containing halogen ions such as
Cl
-
and I
-
. The interactions of ice to various types of ion and its ability preferentially incorporate
certain types of ions depending of ionic species and/or size is a unique characteristic observed by
EIS measurements which can be employed for purposed of ice detection and ionic presence.

48

Figure 3.4.1a: Admittance plots of pure ice using
parallel plate probes.

Figure 3.4.1c: Relative Permittivity plots of pure ice
using parallel plate probes.

Figure 3.4.1e: Loss Tangent plots of pure ice using
parallel plate probes.

Figure 3.4.1b: Nyquist plots of pure ice using
parallel plate probes.

Figure 3.4.1d: Conductivity plots of pure ice using
parallel plate probes.

49

Figure 3.4.2a: Admittance plots of pure ice using
EPC point probes.

Figure 3.4.2c: Relative Permittivity plots of pure ice
using EPC point probes.

Figure 3.4.2e: Loss Tangent plots of pure ice using
EPC point plate probes.

Figure 3.4.2b: Nyquist plots of pure ice using EPC
point plate probes.

Figure 3.4.2d: Conductivity plots of pure ice using
EPC point probes.

50

Figure 3.4.3a: Admittance plots of ice doped with 1
μM KCl using EPC point probes.

Figure 3.4.3c: Relative Permittivity plots of ice
doped with 1 μM KCl using EPC point probes.

Figure 3.4.3e: Loss Tangent plots of ice doped with
1 μM KCl using EPC point plate probes.

Figure 3.4.3b: Nyquist plots of ice doped with 1 μM
KCl using EPC point plate probes.

Figure 3.4.3d: Conductivity plots of ice doped with
1 μM KCl using EPC point probes.

51

Figure 3.4.4a: Admittance plots of ice doped with
100 μM KCl using EPC point probes.

Figure 3.4.4c: Relative Permittivity plots of ice
doped with 100 μM KCl using EPC point probes.

Figure 3.4.4e: Loss Tangent plots of ice doped with
100 μM KCl using EPC point plate probes.

Figure 3.4.4b: Nyquist plots of ice doped with 100
μM KCl using EPC point plate probes.

Figure 3.4.4d: Conductivity plots of ice doped with
100 μM KCl using EPC point probes.

52

Figure 3.4.5a: Admittance plots of ice doped with 1
mM KCl using EPC point probes.

Figure 3.4.5c: Relative Permittivity plots of ice
doped with 1 mM KCl using EPC point probes.

Figure 3.4.5e: Loss Tangent plots of ice doped with
1 mM KCl using EPC point plate probes.

Figure 3.4.5b: Nyquist plots of ice doped with 1 mM
KCl using EPC point plate probes.

Figure 3.4.5d: Conductivity plots of ice doped with 1
mM KCl using EPC point probes.

53

Figure 3.4.6a: Admittance plots of ice doped with 10
mM KCl using EPC point probes.

Figure 3.4.6c: Relative Permittivity plots of ice
doped with 10 mM KCl using EPC point probes.

Figure 3.4.6e: Loss Tangent plots of ice doped with
10 mM KCl using EPC point plate probes.

Figure 3.4.6b: Nyquist plots of ice doped with 10
mM KCl using EPC point plate probes

Figure 3.4.6d: Conductivity plots of ice doped with
10 mM KCl using EPC point probes.

54

Figure 3.4.7a: Admittance plots of ice doped with
100 mM KCl using EPC point probes.

Figure 3.4.7c: Relative Permittivity plots of ice
doped with 100 mM KCl using EPC point probes.

Figure 3.4.7e: Loss Tangent plots of ice doped with
100 mM KCl using EPC point plate probes.

Figure 3.4.7b: Nyquist plots of ice doped with 100
mM KCl using EPC point plate probes.

Figure 3.4.7d: Conductivity plots of ice doped with
100 mM KCl using EPC point probes.

55

Figure 3.4.8a: Admittance plots of ice doped with 1
M KCl using EPC point probes.

Figure 3.4.8c: Relative Permittivity plots of ice
doped with 1 M KCl using EPC point probes.

Figure 3.4.8e: Loss Tangent plots of ice doped with
1 M KCl using EPC point plate probes.

Figure 3.4.8b: Nyquist plots of ice doped with 1 M
KCl using EPC point plate probes.

Figure 3.4.8d: Conductivity plots of ice doped with 1
M KCl using EPC point probes.

56

Figure 3.4.9a: Admittance plots of ice doped with 1
mM KOH using EPC point probes.

Figure 3.4.9c: Relative Permittivity plots of ice
doped with 1 mM KOH using EPC point probes.

Figure 3.4.9e: Loss Tangent plots of ice doped with
1 mM KOH using EPC point probes.

Figure 3.4.9b: Nyquist plots of ice doped with 1 mM
KOH using EPC point probes.

Figure 3.4.9d: Conductivity plots of ice doped with 1
mM KOH using EPC point probes.

57

Figure 3.4.10a: Admittance plots of ice doped with
10 mM KOH using EPC point probes.

Figure 3.4.10c: Relative Permittivity plots of ice
doped with 10 mM KOH using EPC point probes.

Figure 3.4.10e: Loss Tangent plots of ice doped with
10 mM KOH using EPC point probes.

Figure 3.4.10b: Nyquist plots of ice doped with 10
mM KOH using EPC point probes.

Figure 3.4.10d: Conductivity plots of ice doped with
10 mM KOH using EPC point probes.

58

Figure 3.4.11a: Admittance plots of ice doped with
100 mM KOH using EPC point probes.

Figure 3.4.11c: Relative Permittivity plots of ice
doped with 100 mM KOH using EPC point probes.

Figure 3.4.11e: Loss Tangent plots of ice doped with
100 mM KOH using EPC point probes.

Figure 3.4.11b: Nyquist plots of ice doped with 100
mM KOH using EPC point probes

Figure 3.4.11d: Conductivity plots of ice doped with
100 mM KOH using EPC point probes.

59

Figure 3.4.12a: Admittance plots of ice doped with 1
mM CaCl
2
using EPC point probes.

Figure 3.4.12c: Relative Permittivity plots of ice
doped with 1 mM CaCl
2
using EPC point probes.

Figure 3.4.12c: Loss Tangent plots of ice doped with
1 mM CaCl
2
using EPC point probes.

Figure 3.4.12b: Nyquist plots of ice doped with 1
mM CaCl
2
using EPC point probes.

Figure 3.4.12d: Condutivity plots of ice doped with 1
mM CaCl
2
using EPC point probes.

60

Figure 3.4.13a: Admittance plots of ice doped with
10 mM CaCl
2
using EPC point probes.

Figure 3.4.13c: Relative Permittivity plots of ice
doped with 10 mM CaCl
2
using EPC point probes.

Figure 3.4.13c: Loss Tangent plots of ice doped with
10 mM CaCl
2
using EPC point probes.

Figure 3.4.13b: Nyquist plots of ice doped with 10
mM CaCl
2
using EPC point probes.

Figure 3.4.13d: Conducitivity plots of ice doped with
10 mM CaCl
2
using EPC point probes.

61

Figure 3.4.14a: Admittance plots of ice doped with
100 mM CaCl
2
using EPC point probes.

Figure 3.4.14c: Relative Permittivity plots of ice
doped with 100 mM CaCl
2
using EPC point probes.

Figure 3.4.14e: Loss Tangent plots of ice doped with
100 mM CaCl
2
using EPC point probes.

Figure 3.4.14b: Nyquist plots of ice doped with 100
mM CaCl
2
using EPC point probes.

Figure 3.4.14d: Conductivity plots of ice doped with
100 mM CaCl
2
using EPC point probes.

62

Figure 3.4.15a: Admittance plots of ice doped with
100 mM KI

using EPC point probes.

Figure 3.4.15c: Relative Permittivity plots of ice
doped with 100 mM KI

using EPC point probes.

Figure 3.4.15e: Loss Tangent plots of ice doped with
100 mM KI

using EPC point probes.

Figure 3.4.15b: Nyquist plots of ice doped with 100
mM KI

using EPC point probes.

Figure 3.4.15c: Conductivity plots of ice doped with
100 mM KI

using EPC point probes.

63

Figure 3.4.16a: Admittance plots of ice doped with
100 mM MgSO
4
using EPC point probes.

Figure 3.4.16c: Relative Permittivity plots of ice
doped with 100 mM MgSO
4
using EPC point probes.

Figure 3.4.16e: Loss Tangent plots of ice doped with
100 mM MgSO
4
using EPC point probes.

Figure 3.4.16b: Nyquist plots of ice doped with 100
mM MgSO
4
using EPC point probes.

Figure 3.4.16d: Conductivity plots of ice doped with
100 mM MgSO
4
using EPC point probes.

64

Figure 3.4.17a: Admittance plots of ice doped with
100 mM NaCl

using EPC point probes.

Figure 3.4.17c: Relative permittivity plots of ice
doped with 100 mM NaCl

using EPC point probes.

Figure 3.4.17e: Loss Tangent plots of ice doped with
100 mM NaCl

using EPC point probes.

Figure 3.4.17b: Nyquist plots of ice doped with 100
mM NaCl

using EPC point probes.

Figure 3.4.17d: Conductivity plots of ice doped with
100 mM NaCl

using EPC point probes.

65

3.5 Equivalent Circuit Modeling on Ice
Upon the application of the electric field to any insulating material, the resulting polarization
can be categorized into two parts in accordance to the time constant of the response (Barsoukov
and ross Macdonald, 2005, Chpt. 2, pg. 30):
1. An almost instantaneous polarization due to displacement of the electrons with respect to
the nuclei and small distortions such as ionic vibrations of the molecules under the
restoring forces. Polarization of this type defines the high-frequency dielectric constant,
ε
∞
.
2. A time-dependent polarization due to the orientation of dipoles in the electric field. This
case of polarization occurs if the field remains in place for an infinitely long time, and
defines the static dielectric constant, ε
s
.
In the case of ice, however, another process of steady current flow in accordance with Ohm’s
law is also known to occur under suitable electrode configurations. However, this observed
current is attributed to the flow of protons since there is no detectable electronic conduction.
Such a unique conduction process has features in common with electronic conduction in
semiconductors in which ice is aptly referred to as “protonic semiconductors,” though the
protons do not share the quantum mechanical properties that are an essential feature of the band
theory of electrons in solids (Petrenko and Whitworth 1999, pp. 60-61). All three equivalent
circuit models described in the following adequately capture such processes from experimental
measurements on ice.

3.5.1 The Debye Equivalent Circuit
The basis for equivalent circuit modeling in ice is described by the Debye theory in
section 3.1.1. This model is shown in Figure 3.5.1.1., and represents the ideal behavior of any
dielectric materials exhibiting a single time constant.

66

Figure 3.5.1.1 – The Debye Equivalent Circuit Model

The model comprise of only three circuit components, with capacitor 1 (C
1
) in parallel with
resistor 2 (R
2
), and capacitor (C
2
). Mathematical derivation details of the Debye circuit are in
Appendix C.1. The first principles derivation of the Debye equations is in Appendix B. From
admittance measurements, the complex admittance form of the model is as described by equation
3.5.1.1.









 


2
2
2
2
2
2
1
2
2
2
2
2
2
2 2
2
1 1 C R
C
C i
C R
C R
Y




3.5.1.1

The circuit components are re-written as follows:

 
2 2
0 2
0 1
C R
Geo C
Geo C
D
s

 




  
 
3.5.1.2
where Geo = A/L is the geometric factor for the probes configuration and
0
 has a value of
8.85·10
-14
F/cm. Therefore, the equivalent complex relative permittivity form of equation 3.5.1.1
is
 
2 2 2 2
*
1 1
D
s
D
s
j
 
  
 
 
 





 
 

3.5.1.3
From equation 3.5.1.3, real part of the complex relative permittivity equation is the dielectric
constant while the imaginary component is the dielectric loss which is equivalent to the
conductivity (σ) by the following equation 3.5.1.4.

67

Geo
o
"
    3.5.1.4

The Debye relaxation can also be described by the loss tangent, which is the ratio of the
imaginary over the real components of the complex relative permittivity as shown by equation
3.5.1.5.
2 2 '
"
) (
) (
D s
s
Tan
   
  







  3.5.1.5

The model frequency dispersion profiles of the Nyquist, relative permittivity,
conductivity, and loss tangent plots from the Debye model are shown in Figures 3.5.1.1 to
3.5.1.4. For ice, the values of ε
∞
and ε
s
are fixed with only the damping factor of R
2
changing as
a function of both ice structure and temperature, and therefore, controls the frequency response
of the Debye model. The plots of the conductivity and dielectric constant exhibit unique
frequency dispersion which can be used to identify presence of ice and characterize the
microscopic properties of its structure. The magnitude of the loss tangent, however, is
independent of R
2
, with a peak value at 2.5. Overall, the characteristic time constant remains the
most critical factor to understanding the microscopic properties of ice using the Debye model.

68

Figure 3.5.1.1: Cole – Cole Plot for Admittance for
the Debye Model. ε
∞
= 3, ε
s
= 80. τ = 0.0005 sec
(Est 1), τ = 0.0001 sec (Est 2), τ = 0.00005 sec (Est
3).

Figure 3.5.1.3: Conductivity Plot for the Debye
Model. ε
∞
= 3, ε
s
= 80. τ = 0.0005 sec (Est 1), τ =
0.0001 sec (Est 2), τ = 0.00005 sec (Est 3).

Figure 3.5.1.2: Relative Permittivity Plot for the
Debye Model. ε
∞
= 3, ε
s
= 80. τ = 0.0005 sec (Est 1),
τ = 0.0001 sec (Est 2), τ = 0.00005 sec (Est 3).

Figure 3.5.1.3: Loss Tangent Plot for the Debye
Model. ε
∞
= 3, ε
s
= 80. τ = 0.0005 sec (Est 1), τ =
0.0001 sec (Est 2), τ = 0.00005 sec (Est 3).

69

3.5.2 The Debye Equivalent Circuit with Leakage
To account for the proton conduction, the Debye circuit must be modified with the
addition of R
o
circuit component as shown in Figure 3.5.2.1. The R
o
circuit component is often
referred to as the leakage current. Mathematical derivation details of the Debye circuit with
leakage are in Appendix C.2.

Figure 3.5.2.1 – Debye Equivalent Circuit Model with Leakage Current.

The admittance transfer function of the equivalent circuit with leakage is








 

 
2 2
2
2
2
2
1 2 2
2
2
2
2 2
2 2
0
1 1
1




C R
C
C j
C R
C R
R
Y 3.5.2.1

which is virtually identical to the initial Debye equation 3.5.1.1. The only exception is the
additional R
0
term which can be rewritten in terms of conductivity by following

Geo R
DC
 / 1
0
 3.5.2.2

where σ
DC
is represents the low frequency conductivity often referred to as the DC conductivity.
Again, substituting for the same C
1
, C
2
, and τ, along with σ
DC
yields

 














 
 

0
2 2 2 2
*
1 1 

 
  
 
 
 
DC
D
s
D
s
j 3.5.2.3

70

which is the complex relative permittivity form of the equation 3.5.2.1. The loss tangent of the
Debye circuit with leakage is:

   
2
2 2
1
0
2 2
1 2
2
1
/ / 1
'
"
) tan(
C C
R R C
D
 
 
 
  
   


 3.5.2.4

The model frequency dispersion profiles of the Nyquist, relative permittivity,
conductivity, and loss tangent plots from the Debye model with leakage are shown in Figures
3.5.2.1 to 3.5.2.4. This model is very similar to the Debye model, with strong deviation only in
the limiting case at low frequencies where the R
o
circuit component has a unique response for
both the conductivity and the loss tangent profiles. Again, these characteristic responses are
presumably associated with proton conduction near the electrode interface. However, many
physical or chemical processes can also realistically account for the leakage, including surface
heating/melting and even electrode polarization.

71

Figure 3.5.2.1: Cole – Cole Plot for Admittance for
the Debye Model with Leakage. ε
∞
= 3, ε
s
= 80. τ =
0.0005 sec (Est 1), τ = 0.0001 sec (Est 2), τ =
0.00005 sec (Est 3). σ
DC
= 1E-12 S/cm (Est 1), σ
DC
=
1E-10 S/cm (Est 2), σ
DC
= 1E-9 S/cm (Est 2).

Figure 3.5.2.3: Conductivity Plot for the Debye
Model with Leakage. ε
∞
= 3, ε
s
= 80. τ = 0.0005 sec
(Est 1), τ = 0.0001 sec (Est 2), τ = 0.00005 sec (Est
3). σ
DC
= 1E-12 S/cm (Est 1), σ
DC
= 1E-10 S/cm (Est
2), σ
DC
= 1E-9 S/cm (Est 3).

Figure 3.5.2.2: Relative Permittivity Plot for the
Debye Model with Leakage. ε
∞
= 3, ε
s
= 80. τ =
0.0005 sec (Est 1), τ = 0.0001 sec (Est 2), τ =
0.00005 sec (Est 3). σ
DC
= 1E-12 S/cm (Est 1), σ
DC
=
1E-10 S/cm (Est 2), σ
DC
= 1E-9 S/cm (Est 3).

Figure 3.5.2.4: Loss Tangent Plot for the Debye
Model with Leakage. ε
∞
= 3, ε
s
= 80. τ = 0.0005 sec
(Est 1), τ = 0.0001 sec (Est 2), τ = 0.00005 sec (Est
3). σ
DC
= 1E-12 S/cm (Est 1), σ
DC
= 1E-10 S/cm (Est
2), σ
DC
= 1E-9 S/cm (Est 3).

72

3.5.3 The Debye Equivalent Circuit with Proton Conduction and Resonance
Simply stated, the Debye circuit for relaxation essentially describes the mechanical
relaxation from two possible configurations of ice and the characteristic time constant between
them. However, such a transition from one configuration to another due to macroscopic elastic
strain from the applied electric field typically has an associated resonance effect; such viscous
damping forces exhibits a maxima at a particular resonance frequency (Fletcher, pp 174). Like
conventional parallel plate dielectric measurements, our unique probe configuration does detect
the resonance. Such resonance should accompany the slow transition in configuration states in
ice by way of energy dissipation in dielectric spectroscopy measurements. To account for the
such resonance effect, we added an inductor component, L
2
, in series with the R
2
and C
2
circuit
components to the Debye equivalent circuit model as shown in Figure 3.5.3.1. Mathematical
derivation details of the Debye circuit with resonance are in Appendix C.3.

Figure 3.5.3.1 – Debye Equivalent Circuit Model and Resonance.

The corresponding admittance transfer function of the equivalent circuit with resonance is

 
 
 








 

 
 
 
2 2
2
2
0
2
2
0
2
2
1
2 2
2
2
0
2
2
2 2
0 / 1
/ 1
/ 1
/ 1
D D
C
C j
R
R
Y
   
 

   
 
3.5.3.1

where τ is the same as the previous circuits, the ω
o
is the new parameter representing the natural
resonance frequency of ice during relaxation as described by equation 3.5.3.2.

73

2 2 0
/ 1 C L  
3.5.3.2

Therefore, the Debye time constant at resonance angular frequency (τ
R
) is defined as

o
R


1

3.5.3.3

Substituting in for ω
o
, the complex relative permittivity form of the model is Equation 3.5.3.3.

   
 
 
 









 


 
 
 
 

0
2 2
2
2
0
2 2 2
2
2
0
2
2
0
2
*
/ 1 / 1
/ 1


   
  
   
   
 
DC
D
s
D
s
j 3.5.3.4

Similarly, the loss tangent equation describing energy dissipation is 3.5.3.6.

   
     
2
0
2
2
2 2
2
2
0
2
1
2
2 2
2 2
2
2
0
2
1
/ 1 / 1
/ 1
1
'
"
) tan(
       
 
   



   
  
 
C C
R R
D
D
D
3.5.3.5

The model frequency dispersion profiles of the Nyquist, relative permittivity, conductivity, and
loss tangent plots from the Debye model with leakage are shown in Figures 3.5.3.1 to 3.5.3.4.
The Nyquist plot with the addition of the inductance circuit component exhibits an inductance
loop the low frequency range. The conductivity curves sag at high frequencies and this is a
manifestation of resonance behavior as well. The relative permittivity curves exhibits the
classical Debye behavior with exception of the dip at the resonance frequency. In this region, the
dielectric constant for high values of inductance has a value less then unity, which is not
physically possible. That is, the dielectric constant cannot be less than unity since the dielectric
constant for a vacuum and very close to the value for air. Again, this behavior is due to
molecular resonance and requires an RLC circuit to model this behavior.
The tangent curves with the resonance effect are also very unique features as compared to
the classical Debye circuit. At low frequencies where 1
2 2

D
  , the log tan θ vs log f has a
slope of +1, and equates to Equation 3.5.3.6.

74

 
D
C C
D
C C
C
D

 

 











1 2
2 2
2 1
2
1
tan 3.5.3.6

At high frequencies where 1
2 2

D
  , the log tan θ vs log f plot (Figure 3.5.3.4) approaches a
slope of -3, and equates to Equation 3.5.3.7.

 
3
1
2
2
2
1
2 2 tan


 
C L
R


3.5.3.7

Also apparent in the loss tangent curves are the increase in the peak height over the 2.5 value
calculated by the Debye circuit without the presence resonance. The loss tangent curves
intersect at tan θ = 1. This is the same behavior observed for the phase angle where the curves
pass through 45°. Note that the inductance influences only the high frequency side of the loss
tangent. The slope of the high frequency side of the loss tangent increases from a slope of -1 for
L
2
= 0 to a slope approaching -3 for large values of L
2
.

75

Figure 3.5.3.1: Cole – Cole Plot for Admittance for
the Debye Model with Leakage. ε
∞
= 3, ε
s
= 80, τ =
0.00055 sec, σ
DC
= 1E-20 S/cm. ω
o
= 8.53E-3 sec
-1

(Est 1), ω
o
= 9.53E-3 sec
-1
(Est 2), ω
o
= 1.29E-4 sec
-1

(Est 3).

Figure 3.5.3.3: Conductivity Plot for Admittance for
the Debye Model with Leakage. ε
∞
= 3, ε
s
= 80, τ =
0.00055 sec, σ
DC
= 1E-20 S/cm. ω
o
= 8.53E-3 sec
-1

(Est 1), ω
o
= 9.53E-3 sec
-1
(Est 2), ω
o
= 1.29E-4 sec
-1

(Est 3).

Figure 3.5.3.2: Relative Permittivity Plot for
Admittance for the Debye Model with Leakage. ε
∞
=
3, ε
s
= 80, τ = 0.00055 sec, σ
DC
= 1E-20 S/cm. ω
o
=
8.53E-3 sec
-1
(Est 1), ω
o
= 9.53E-3 sec
-1
(Est 2), ω
o
=
1.29E-4 sec
-1
(Est 3).

Figure 3.5.3.4: Loss Tangent Plot for Admittance for
the Debye Model with Leakage. ε
∞
= 3, ε
s
= 80, τ =
0.00055 sec, σ
DC
= 1E-20 S/cm. ω
o
= 8.53E-3 sec
-1

(Est 1), ω
o
= 9.53E-3 sec
-1
(Est 2), ω
o
= 1.29E-4 sec
-1

(Est 3).

76

3.5.4 Equivalent Circuit Modeling Results
A major challenge in equivalent circuits modeling on ice using complex non-linear least
square (CNLS) techniques is attributable to poor data quality due to instrument limitation
especially at low temperatures. To maximize model accuracy, equivalent circuit curve fitting
protocols were conducted using a combination of both ZSimWin (version 3.10) and algorithms
developed from Matlab from Mathworks. ZSimWin, from Princeton Applied Research, Inc.,
performs complex non-linear least square fitting of impedance measurement using a Simplex
algorithm for initial component estimation. Subsequently, the fitting results from ZSimWin were
the initial inputs for the matlab fitting algorithm by constraining the equilibrium polarization
circuit components of C
1
and C
2
which correspond to valid range of magnitudes of ε
∞

(approximately 3) and ε
s
(80-100) values. The circuit component R
o
representing ice proton
conduction is extracted from the high frequency limit, leaving only two circuit components, R
2

and L
2
, as the only unknowns for fitting optimization. Thus, the matlab fitting algorithm is based
on data regression on the individual complex real and imaginary components of the raw
admittance data solving R
2
and L
2
to minimize the error based on





2
2
2
) (
) (
Ymean Yact
Ymean Yest
R 3.5.4.1
where R
2
is the square of the residual error, Yest model admittance, Yact is the measured
admittance, and Ymean is the average value of the admittance.
Since the ideal Debye equivalent circuit does not account for current leakage attributable
to the interfacial liquid layer and/or proton conduction in the low frequencies, the measurement
results from parallel plate electrode configuration is modeled using the modified Debye
equivalent circuit with leakage as shown in Figure 3.5.2.1. A direct comparison of model results
with actual measurements on Nyquist, loss tangent, conductivity, and relative permittivity plots
at ice temperatures ranging from -15
o
C to -65
o
C are shown in Figures 3.5.4.1a-d, respectively.
The Nyquist plot exhibits not inductance loop indicating resonance, and the single Debye
relaxation as defined by radius of the semi-circle decreases with decreasing temperature. As
expected, the addition of the leakage resistance R
o
shifts the intercept to the Y
r
–axis greater than

77

zero at low frequency limit as a function of liquid water at the probe interface. From the loss
tangent spectra, the loss tangent peak value is in excellent agreement at 3 and is relatively
constant at all ice temperatures. However, the peak loss tangent is slightly higher than the
originally predicted value of 2.5. At higher temperatures such as 15
o
C, a second relaxation is
observe and is also represented accurately by the equivalent circuit. The modeled conductivity
spectra exhibited a constant region at the high frequency limit which is classically interpreted as
the region dominated by the thermodynamic Bjerrum protonic defect concentration. The ionic
defects dominate region indicated by the convergence of the sloped region at the low frequency
range. Finally, the modeled relative permittivity is in excellent agreement with the values of ε
∞

and ε
s
, along with the frequency corresponding to the Debye relaxation which indicates no
significant deviation over the entire spectra for all ice temperatures.
Measurements performed using the point electrode configurations also exhibited good
agreement with the modeled results. Here, ice measurements exhibit a distinctive resonance
relaxation which is accurately represented by the equivalent circuit with resonance as shown in
Figure 3.5.3.1. Model comparison plots are shown in Figures 3.5.4.2a-d. From the Nyquist plot,
the model results exhibits an inductive loop like the actual measurements. The loss tangent
spectra from the model also accurately agree with the Debye relaxation peaks at different
temperatures similar to the actual measurements. These relaxation peak increases from 5 at -15
o
C to approximately 100 at -45
o
C which are a very unique behavior in the peak loss tangent as
compared to the parallel plate probe measurements. Also unique to the point electrode
configuration, the conductivity in the high frequency region exhibits a slight negative slope as
modeled by the equivalent circuit and drops-off precipitously. The relative permittivity model
exhibits the dip corresponding to the Debye resonance similar to the actual measurements.
However, significant deviation at the onset of the resonance frequency is observed and increases
with increasing temperature. This is most likely the result of Maxwell polarization which also
tends to increase with higher temperatures and is attributable to the presence of liquid water layer
at the interface.
Doped ice is also accurately modeled with the Debye equivalent circuit with resonance.
A representative comparison of ice doped with 100 mM KCl is shown in Figures 3.5.4.3a-d. The
inductance loop exhibits a longer tail as a consequence of the increased in Debye relaxation

78

frequency as observed from the Nyquist plots. The loss tangent peak is significantly smaller than
that of pure ice ranging from 4 to 7 from -25
o
C to -65
o
C, respectively. Like pure ice, the
modeled loss tangent of 100 mM KCl doped ice also exhibits significant deviation at the
resonance frequency. Though the loss tangent peak is accurately predicted to increase with
decrease in temperature, the modeled results are lower than that of the actual measurements. The
modeled conductivity of doped ice is similar to pure ice with significant deviation observed at
the high frequencies as well as the low frequency deviation as observed from the relative
permittivity spectra due to Maxwell polarization. These large deviations are mainly attributable
to the L
2
circuit component which tends to be sensitive to both temperature as well as Maxwell
polarization and poses a significant challenge in maximizing model convergence.

79

Figure 3.5.4.1a: Nyquist plot comparing actuals to
equivalent circuit modeling results using the Debye
equivalent circuit with leakage on pure ice. Plate
electrode configuration.

Figure 3.5.4.1c: Conductivity spectra comparing
actuals to equivalent circuit modeling results using
the Debye equivalent circuit with leakage on pure ice.
Plate electrode configuration.

Figure 3.5.4.1b: Loss tangent spectra comparing
actuals to equivalent circuit modeling results using
the Debye equivalent circuit with leakage on pure ice.
Plate electrode configuration.

Figure 3.5.4.1d: Relative permittivity spectra
comparing actuals to equivalent circuit modeling
results using the Debye equivalent circuit with
leakage on pure ice. Plate electrode configuration.

80

Figure 3.5.4.2a: Nyquist plot comparing actuals to
equivalent circuit modeling results using the Debye
equivalent circuit with leakage on pure ice. Point
electrode configuration.

Figure 3.5.4.2c: Conductivity spectra comparing
actuals to equivalent circuit modeling results using
the Debye equivalent circuit with leakage on pure ice.
Point electrode configuration.

Figure 3.5.4.2b: Loss tangent spectra comparing
actuals to equivalent circuit modeling results using
the Debye equivalent circuit with leakage on pure ice.
Point electrode configuration.

Figure 3.5.4.2d: Relative permittivity spectra
comparing actuals to equivalent circuit modeling
results using the Debye equivalent circuit with
leakage on pure ice. Point electrode configuration.

81

Figure 3.5.4.3a: Nyquist plot comparing actuals to
equivalent circuit modeling results using the Debye
equivalent circuit with leakage on pure ice. Point
electrode configuration.

Figure 3.5.4.3c: Conductivity spectra comparing
actuals to equivalent circuit modeling results using
the Debye equivalent circuit with leakage on pure ice.
Point electrode configuration.

Figure 3.5.4.3b: Loss tangent spectra comparing
actuals to equivalent circuit modeling results using
the Debye equivalent circuit with leakage on pure ice.
Point electrode configuration.

Figure 3.5.4.3d: Relative permittivity spectra
comparing actuals to equivalent circuit modeling
results using the Debye equivalent circuit with
leakage on pure ice. Point electrode configuration.

82

The summary of the fitting parameters on the circuit components of the Debye equivalent
circuit with resonance for modeling pure and 100 mM doped ice measurements are listed in
Table 3.5.4.1 for all measured temperatures. No model results are available for measurements
exhibiting liquid-like behavior. Electrolytes containing 100 mM CaCl
2
exhibited ice signatures
at temperatures below -55
o
C. At -15
o
C, only 100 mM KCl electrolytes exhibited ice
characteristics. The resistor component R
o
corresponding to leakage current is consistent for all
dopant types in the range of 10
7
to 10
10
Ω. As suggested earlier, R
o
accounts for proton
conduction from liquid-like behavior at the interface. Thus, it is expected with ice doped with
KCl at -15
o
C should exhibit the lowest value while KI doped ice exhibits the highest value at
approximately 10
10
Ω. Also as expected, the values corresponding to capacitive components C
1

and C
2
is relatively constant at approximately 3x10
-13
F and 5x10
-13
F, which correspond to the
static and infinite frequency dielectric constant of ε
s
and

ε
∞
, respectively. Ice measurements
doped with CaCl
2
and KOH, however, resulted in significantly higher C
1
values in the order of
5x10
-11
F, and is attributable to Maxwell polarization and the presence of water layer at the
interface. The values of the circuit component R
2
is the resistance in the bulk and varies
substantially depending on dopant type and temperature. The lowest value of R
2
is exhibited by
ice doped with KOH at -45
o
C with a value in the order of 10
5
Ω. Pure ice and ice doped with
MgSO
4
exhibited the highest R
2
values at 10
9
Ω, and is consistent with the low concentration of
the larger Mg
+2
and SO
4
-2

ions in the lattice. Substantial variability is also observed in L
2
values.
The resonance is highly dependent on temperature, corresponding to the values of L
2
increasing
with decrease in temperature. Ice doped with halogen ions as well as OH
-
ions tends to exhibit
the lowest L
2
values, while pure ice exhibits exhibit the highest value. The highest L
2
values are
at -65
o
C for pure ice and MgSO
4
doped ice at approximately 10
7
H. Overall, the circuit
components are in the range of expected values representing various physical processes to be
addressed in the following discussions.

83

Table 3.5.4.1: Summary of circuit component fitting results using Debye circuit with resonance
for 100 mM doped ice.

Temperature (
o
C)
Dopant R
0
C
1
R
2
C
2
L
2
-15 None (Pure Water) 3.30E+08 2.60E-13 2.20E+07 4.50E-12 9.00E+01
KCl 1.00E+07 2.60E-13 1.20E+06 4.50E-12 6.00E-01
KOH -- -- -- -- --
*CaCl
2
-- -- -- -- --
KI -- -- -- -- --
MgSO
4
-- -- -- -- --
NaCl -- -- -- -- --
-25 None (Pure Water) 4.20E+08 2.60E-13 6.80E+07 4.50E-12 8.50E+02
KCl 8.00E+08 3.40E-13 2.65E+06 5.60E-12 2.00E+00
KOH -- -- -- -- --
*CaCl
2
-- -- -- -- --
KI -- -- -- -- --
MgSO
4
1.43E+08 2.70E-13 7.00E+07 5.00E-12 1.35E+03
NaCl -- -- -- -- --
-35 None (Pure Water) 4.70E+08 2.60E-13 2.20E+08 4.50E-12 1.10E+04
KCl 1.00E+09 3.00E-13 3.90E+06 5.60E-12 3.50E+00
KOH -- -- -- -- --
*CaCl
2
-- -- -- -- --
KI 1.34E+09 2.56E-13 2.86E+06 6.00E-12 1.40E+00
MgSO
4
1.91E+08 2.70E-13 2.60E+08 5.00E-12 2.80E+04
NaCl 1.73E+09 2.17E-13 4.01E+06 5.00E-12 2.70E+00
-45 None (Pure Water) 6.20E+08 2.60E-13 1.00E+09 4.50E-12 3.70E+05
KCl 4.00E+09 3.00E-13 7.04E+06 5.84E-12 1.20E+01
KOH 3.67E+05 3.69E-13 8.81E+05 3.00E-11 1.00E-01
*CaCl
2
-- -- -- -- --
KI 4.97E+09 2.58E-13 3.70E+06 5.50E-12 2.50E+00
MgSO
4
2.72E+08 2.30E-13 1.00E+09 6.00E-12 5.10E+05
NaCl 2.20E+09 2.24E-13 6.66E+06 5.00E-12 8.00E+00
-55 None (Pure Water) 9.58E+08 2.44E-13 2.55E+09 4.00E-12 4.50E+06
KCl 4.00E+09 2.70E-13 1.21E+07 6.15E-12 3.00E+01
KOH 2.80E+06 3.62E-13 2.27E+06 1.50E-11 1.00E+00
*CaCl
2
1.36E+08 2.72E-13 6.71E+07 2.00E-11 1.40E+03
KI 9.73E+09 2.54E-13 4.96E+06 5.50E-12 4.20E+00
MgSO
4
5.58E+08 2.30E-13 1.50E+09 6.00E-12 1.10E+06
NaCl 3.70E+09 2.26E-13 1.30E+07 5.00E-12 3.10E+01
-65 None (Pure Water) 1.86E+09 2.43E-13 4.49E+09 4.00E-12 8.50E+06
KCl 8.00E+09 2.70E-13 2.50E+07 6.15E-12 1.35E+02
KOH 1.90E+07 2.53E-13 6.52E+06 9.00E-12 6.00E+00
*CaCl
2
2.36E+08 2.72E-13 1.01E+08 2.50E-11 2.85E+03
KI 4.11E+10 2.61E-13 7.24E+06 5.50E-12 9.00E+00
MgSO
4
9.78E+08 2.30E-13 3.30E+09 2.00E-12 5.21E+06
NaCl 4.30E+09 2.24E-13 2.69E+07 5.00E-12 1.35E+02
Units: R (=) Ohms, C (=) Farads, L (=) Henrys
*CaCl
2
circuit fitting analysis was conducted with measurements at temperature range between -55 to -70
o
C.

84

3.6 Results and Discussion

The crystalline nature of ice is such that it behaves like a semiconductor, exhibiting
properties of both dielectric and conductive material. Therefore, the inherent electrical
properties as measured by EIS technique are unique in ice in terms of polarization, relaxation,
and conductive properties. Combined with the equivalent circuit modeling which reasonably
describes the measurement results, direct understanding of the electrical behavior emanating
from the molecular level can be quantitatively achieved from measurement parameters such as
relaxation times, resonance times, and bulk conductivity of the ice crystal.

3.6.1. Debye Relaxation Times
Ice exhibits a single relaxation process aptly name Debye relaxation time constant which
can be directly calculated from fitted values of the circuit components achieved through
modeling protocol described in the earlier section. Figure 3.6.1.1 compares the relaxation times
of the ideal (plate electrodes) and vs. point electrode (fieldable) configurations over the same
temperature range. At temperatures up to -20
o
C, the time constants from both point and plate
electrode configuration were in the range of 10
-4
sec. This is in excellent agreement with
numerous past studies conducted at similar conditions.
46,47
The point electrodes deviation from
the ideal cases is only significant at low temperatures below -40
o
C with deviation greater than
20% and can be attributed to the both instrument limit as well as electrode design. With values
of relaxation times, the activation energy (E
a
) for Debye relaxation can be ascertained using the
form of the Arrhenius equation







 
kT
E
A
a
D
exp  3.6.1.1
where A is the pre-exponential coefficient in seconds, k is the Boltzmann constant, and T is the
temperature in kelvin. The –ln(τ
D
) vs 1/T plot comparing the two electrode configurations is
shown in Figure 3.6.1.2 in which the value of E
a
is extracted from their respective slopes. The
activation energy for the plate and point electrodes are determined at 0.529 eV and 0.579 eV,
respectively. These values for pure ice are also in good agreement with to values reported in
earlier studies (Haas, 1962) and represent the intrinsic behavior of pure ice.

85

Figure 3.6.1.1: Comparison of Debye relaxation time constants from measurements using plate
electrodes vs point electrodes on pure ice.

Figure 3.6.1.2: Arrhenius plots of Debye relaxation time constants from plate electrodes vs point
electrodes on pure ice. Activation energies for Debye relaxation: parallel plate electrode
configuration = 0.529 eV, point electrode configuration = 0.579 eV.

86

The summary of Debye relaxation times for KCl doped ice at various electrolyte
concentrations is shown in Figure 3.6.1.3. At a molecular level, the relaxation process in ice
directly corresponds to the electrical property of susceptibility strength. As expected, the Debye
relaxation time in ice is extremely sensitive to dopants such as K
+
and Cl
-
ions, and generally
increases with increasing concentration and decrease with temperature. The only exception is
the values reported above -40
o
C with 1 mM dopant concentration, where values increased as
much as 1 order of magnitude at ice temperatures near -15
o
C. This suggests the mechanisms for
relaxation in pure ice is different than that observed in doped ice and can be attributable to their
activation energies associated with defect formation and transport. At concentrations of 10 mM
and above, the time constants for the Debye relaxation decreases as expected by as much as 2
orders of magnitude lower. The lowest measured value is at -15
o
C at 5.4 x 10
-6
sec for 100 mM
KCl, which supports the basic understanding of lower relaxation times due to lower crystallinity
and higher defect concentration. The decrease in relaxation times is inversely proportional to the
dopant concentration up to 100 mM. The time constants for 100 mM and 1 M dopant
concentration of K
+
and Cl
-
ions are very similar, and are attributable to the saturation
concentration of ion absorption into the ice lattice. Thus, conventional theory that asserts defect
concentration is increased by ion absorption into the ice lattice is upheld as observed by
decreases in the Debye relaxation times.
From equation 3.6.1.1, the activation energy is determined from the Arrhenius plot at
various KCl dopant concentrations in the ice. Table 3.6.1.1 summarizes the activation energies
(units of eV) of various concentrations KCl doped ice as compared to pure ice. The ice doped
with KCl salt is roughly 50% lower than that of pure ice at approximately 0.30 eV. The
relaxation activation energy is inversely proportional to the KCl dopant concentration up to a
saturation concentration of 100 mM. For 1 mM, a value of 0.20 eV is found as compared to a
value of 0.29 eV at 100 mM and beyond. Lower activation energies suggest lower energy
barriers of defect formation and migration. The observed lower relaxation times at the highest
KCl dopant concentration of 100 mM suggest defect concentration remains the dominating factor
in the Debye relaxation mechanism and its measured relaxation times.

87

Figure 3.6.1.3: Comparison of Debye relaxation time constants of KCl doped ice at various
electrolyte concentrations.

Table 3.6.1.1: Summary of E
a
values from Debye relaxation times (τ
D
) at various dopant
concentrations of KCl.

The Debye relaxation time constants were also compared against various dopant types up
to an ionic concentration of 100 mM. The summary plot is shown in Figure 3.6.1.4. Relative to
the time constants of pure ice, the doped ice from MgSO
4
electrolyte exhibited the most similar
Debye relaxation times over the entire ice temperature range from -15
o
C to -65
o
C. These
similarities suggest that larger ions such SO
4
-2
is relatively insoluble in the ice lattice despite the
high valence number. The lowest time constant values of doped ice are those from KI and KOH
electrolytes. However, it is of note that ice doped with KOH may exhibited higher time
KCl Dopant
Concentration E
a
(eV)
Pure Ice 0.579
1 mM 0.199
10 mM 0.273
100 mM 0.293
1 M 0.292

88

constants than expected as a consequence of model error attributable to the strong presence of
Maxwell polarization at the interface. Therefore, both halogen and OH
-
ions exhibit similar
solubility in ice which can be attributable to the disorder increases to the hydrogen network upon
ion absorption. In addition, ice doped with KCl and NaCl exhibited very similar values of Debye
time constants and highlights the similarities in their freezing mechanisms in the presence of
such ions. The ice doped with CaCl
2
electrolyte, however, exhibited the second highest values of
time constants at freezing temperature below -60
o
C, which suggests that an excess of the same
halogen ions in combination with larger ion does not increase ion absorption into the ice lattice.
The size of cation appears to play a significant role in the Debye relaxation mechanism, and has
a tendency to decrease the relaxation times with larger cations more than smaller anions.
A comparison of activation energies from various dopant types was also calculated as
shown in Table 3.6.1.2. The lowest activation energy is exhibited by the KI electrolyte at 0.154
eV. Similarly, the ice doped with KOH is determined to be 0.164 eV. These results are
consistent with the work of Kawada (1989). The activation energy for KCl and NaCl are similar
at 0.293 eV and 0.273 eV, respectively. The activation energy from CaCl
2
electrolytes is slightly
higher than that of KCl and NaCl at 0.307 eV, which again highlights the higher solubility of the
halogen anion in ice during the freezing process. The activation energy from MgSO
4
salt is
nearly the same as pure ice at 0.522 eV, indicating very little or no ion concentration in the ice
lattice by Mg
+2
or SO
4
-2
ions leading to higher phase separation due to ion-pairing of the MgSO
4

salt in the ice. The activation energies of the Debye relaxation are highly dependent on ion type,
and thus, can provide significant details on relative ion absorption capabilities of certain ionic
species as well as the intrinsic behavior of the ice lattice structure.

89

Figure 3.6.1.4: Comparison of Debye relaxation time constants of various types of ionic dopants
in ice at 100 mM electrolyte concentrations. The relaxation times of CaCl
2
doped ice is
measured over lower temperatures to ensure the presence of ice.

Table 3.6.1.2: Summary of E
a
values from Debye relaxation times (τ
D
) of various dopant types
at 100 mM dopant concentrations.

Dopant Type E
a
(eV)
Pure Water 0.579
KCl 0.293
NaCl 0.273
KI 0.154
MgSO
4
0.522
CaCl
2
0.307
KOH 0.164

90

3.6.2 Debye Resonance Times
Like the relaxation times, the Debye resonances times were also determined from the
equivalent circuit model parameters as defined by Equation 3.5.3.3. This resonance time
constant is unique to the point electrode configuration and is not observed in the parallel plate
configurations employed by previous researchers. As with the Debye relaxation process, the
activation energy for the Debye resonance can be determined by






 
kT
E
A
a
R
exp  3.6.2.1
where E
a
is the activation for the resonance observed in ice. Figure 3.6.2.1 summarizes the
resonance time constants of KCl doped ice from various ionic concentrations. Similar to the
Debye relaxation times, the magnitude of the resonance times also decreases with increasing KCl
concentrations by more than 2 orders of magnitude. The resonance time constants from 100 mM
and 1 M KCl doped ice are also similar over the measured temperature range, which again
indicates a saturation concentration for ionic absorption. At 1 mM dopant concentration, a
crossover with the pure ice is observed at -25
o
C, and is similar that observed with the Debye
relaxation time constants. The higher resonance time constants relative to pure ice can also be
attributable to differences in activation energies associated defect formation. This is supported
by the Arrhenius plot of –ln(τ
R
) vs 1/T as shown in Figure 3.6.2.2, which illustrates the
differences in the slope of pure ice compared to KCl doped ice. As summarize in Table 3.6.2.2,
the activation energy for resonance in pure ice is 0.667 eV; the activation energy for resonance
for KCl doped ice is substantially lower in the range of 0.20 to 0.30 eV, and incrementally
increases with increasing dopant concentration up to the saturation concentration of 100 mM.
These activation energy values are also similar to that observed for the Debye relaxation, and
strongly suggest a strong similarity in the mechanisms which govern their respective processes.

91

Figure 3.6.2.1: Comparison of Debye resonance time constants of KCl doped ice at various
electrolyte concentrations.

Figure 3.6.2.2: Arrhenius plots of Debye resonance time constants of KCl doped ice at various
electrolyte concentrations.

92

Table 3.6.2.1: Summary of E
a
values from Debye resonance times (τ
R
) at various dopant
concentrations of KCl.

The corresponding resonance times from various dopant types are summarized in Figure
3.6.2.3. As expected, the trends are similar to that of Debye relaxation times. KOH doped ices
exhibited the lowest resonance time constants in the order of 10
-6
sec. Ice doped with NaCl and
KCl is similar in the range of 10
-5
sec. However unlike the Debye relaxation times, KI doped
ices exhibited higher resonance time constants than that of KCl and NaCl doped ice. The values
for CaCl
2
doped ice are in the expected range of 10
-4
sec. The pure ice and MgSO
4
doped ice
exhibited the highest values in the range of 10
-4
to 10
-3
sec over the temperature range of -15
o
C
to -65
o
C, respectively. The activation energies for various types of doped ice are summarized in
Table 3.6.2.2. Here, the lowest activation energy for resonance is exhibited by KI doped ice at
0.125 eV. Ice doped with CaCl
2
and NaCl, and KCl is nearly identical suggesting halogenated
ions have the strongest effect on the resonance process. The highest activation energy for
resonance is determined for MgSO
4
doped ice at 0.549 eV. Therefore, both ion type and size
significantly impact lattice structure. The Debye resonance, similar to the Debye relaxation, is
also dependent on ionic type.
KCl Dopant
Concentration E
a
(eV)
Pure Ice 0.667
1 mM 0.197
10 mM 0.270
100 mM 0.254
1 M 0.283

93

Figure 3.6.2.3: Comparison of Debye resonance time constants of various types of ionic dopants
in ice at 100 mM electrolyte concentrations.

Table 3.6.2.2: Summary of E
a
values from Debye resonance times (τ
R
) of various dopant types
at 100 mM dopant concentrations.

Dopant Type E
a
(eV)
Pure Water 0.667
KCl 0.254
NaCl 0.280
KI 0.125
MgSO
4
0.549
CaCl
2
0.280
KOH 0.295

94

3.6.3 Ice Conduction
The most important parameter for understanding the transport properties of ice is the
conductivity. With the measured conductivity, quantitative information on defect concentration
and mobility can be ascertained from the corresponding ionic species. Ionic defects are not
thermally activated and thereby dominate the low frequency range of the conductivity spectra.
The classical basis for conduction by ionic protonic defects is the probabilistic mechanism of
proton hopping from occupied site to vacant adjacent sites (Petrenko and Whitworth 1999, pp.
131). Ionic defects should also have a profound impact on the Maxwell polarization processes
which typically occur adjacent to the electrodes due to the migration and accumulation of such
charge carrying defects from an applied electric field. On the other hand, the bulk behavior of
ice is thermally activated and comes from the conductivity value discounting all build-up of
polarization in the high frequencies, which means the characterization of Bjerrum type defects is
also made in the high frequency region of the measured spectra.
Similar to the Debye relaxation, the conductivity measurements of pure ice from the point
electrode configuration is compared to that of the standard parallel electrodes as shown in Figure
3.6.3.1. The temperature dependence of the conductivity on ice follows a monotonic behavior
with a high of about 10
-7
S/cm at -15
o
C and corresponds to the similar value of liquid water
prior to freezing; the lower conductivity value of about 10
-10
S/cm at the lowest temperature of -
65
o
C confirms the thermal activation of the conductivity most likely due to motion of charge
carriers such as protonic Bjerrum defects. As expected, the conductivity of pure ice conforms to
the Arrhenius equation shown on Figure 3.6.3.2 with a lower slope exhibited by the parallel plate
electrode configuration as compared to the point electrode configuration. Thus, the
corresponding activation energy for the parallel plate electrode configuration is determined to be
significantly lower at 0.456 eV as compared to the activation energy of 0.543 eV for the point
electrode configuration. This difference is also attributable to the orientation of the applied
electric field due to the point electrode configuration.

95

Figure 3.6.3.1: Comparison of conductivity over the ice temperature range of -15
o
C to -65
o
C
from measurements using plate electrodes vs point electrodes on pure ice. Data were extracted at
10 kHz frequency.

Figure 3.6.3.2: Arrhenius plots of conductivity over the ice temperature range of -15
o
C to -65
o
C from measurements using plate electrodes vs point electrodes on pure ice. Activation
energies for conduction: parallel plate electrode configuration = 0.456 eV, point electrode
configuration = 0.543 eV.

96

In the high frequencies absent of polarization effects, the conductivity of ice is dependent
on the ion concentration. The comparative summary of the conductivity of KCl doped ice at
ionic concentrations ranging from 1 mM to 1 M is shown in Figure 3.6.3.3. As expected, the
temperature dependency is the small and increases with temperature. However, the ice
conductivity of even 1 mM of KCl dopant resulted in as much as 1-2 orders of magnitude higher
(~10
-6
S/cm). Increasing KCl dopant concentration does not proportionally increase the ice
conductivity, and again suggest an ion saturation concentration in the ice. The corresponding
activation energy for conductivity is relatively constant with increasing KCl concentration.
Compared to the pure ice value of 0.543 eV for pure ice, the activation energy of KCl doped ice
is substantially lower in the range of 0.350 +/- 0.001 eV. Though this strongly suggests that
charge carries in KCl doped ice are much more mobile than that of pure ice, it is yet to be
determined that the enhancement in observed conductivity is attributable to factors such as
increased ion or defect concentration.
Similarly, the comparison of various dopant types on the conductivity was also
performed as shown in Figure 3.6.3.4. The highest conductivity is exhibited by CaCl
2
and KOH
doped ice with an increase of nearly 3 orders of magnitude (~10
-5
S/cm) over pure ice. Such
high values are beyond that of liquid water and contradicts the theory that defect charge carriers
as the sole mechanism for ice conduction. Alternatively, the only other viable explanation is the
free ion motion in the ice lattice as oppose to immobile ions trapped in ice cages. The
conductivity of ice doped with halogenated salts exhibited similar increases in conductivity of 1-
2 orders of magnitude over pure ice. In these cases, such substantial increases in conductivity
are most likely attributable to increases in either protonic defects or ionic transport. Lastly, the
conductivities of MgSO
4
doped ice are approximately same as pure ice as a consequence of the
ion-pairing during ice nucleation. Therefore, ions in ice can exhibit a variety of behavior
depending on its physical properties, including immobility or trapped in the ice cavities, free
migration similar to ion motion in liquids, and recombine back to its salt causing phase
separation.
The mechanisms involved conductivity of various dopant types can also be further
examined by their individual activation energies. Table 3.6.3.2 summarizes the activation
energies of various dopant species. Larger ions appear to exhibit the lowest activation energies,

97

with KI doped ice exhibiting the lowest value at 0.176 eV. Ice doped with MgSO
4
exhibited the
next lowest activation energy of 0.309 eV. The highest values of activation energy for
conductivity were exhibited by KOH and CaCl
2
doped ice at 0.495 eV and 0.391 eV,
respectively. Such unique behavior can be explained by ion interaction and transport behavior
which play a significant role in the conductivities of both liquids and ice. Therefore, our
investigation strongly suggests ice doped KOH and CaCl
2
must also contain free ion motion
similar to liquid water which is not significantly observed in ice doped with KCl, NaCl, KI, or
MgSO
4
.

98

Figure 3.6.3.3: Comparison of Debye relaxation time constants of KCl doped ice at various
electrolyte concentrations.

Table 3.6.3.1: Summary of activation energies (E
a
) for conductivity at various dopant
concentrations of KCl.

KCl Dopant
Concentration E
a
(eV)
Pure Ice 0.543
1 mM 0.354
10 mM 0.363
100 mM 0.343
1 M 0.368

99

Figure 3.6.3.4: Comparison of Debye relaxation time constants of various types of ionic dopants
in ice at 100 mM electrolyte concentrations.

Table 3.6.3.2: Summary of activation energies (E
a
) values for conductivity of various dopant
types at 100 mM dopant concentrations.

Dopant Type E
a
(eV)
Pure Water 0.543
KCl 0.343
NaCl 0.368
KI 0.176
MgSO
4
0.309
CaCl
2
0.391
KOH 0.495

100

Finally, a method of quantifying the ion concentration in ice is a critical step towards the
understanding of ice transport processes. As described by Equation 3.1.3.1, the bulk
conductivity of doped ice can provide estimations of the ion concentration. From classical
theories, this quantity is more accurately interpreted as the combined ion concentration and the
corresponding defect count in ice since they are indistinguishable by bulk measurements such as
EIS. Our estimation of the ion concentration in ice employs the calibration data of conductivity
at various ionic concentrations in the liquid phase assuming free ion motions and a chemical
activity of unity at dilute liquid concentrations. In particular the case of KCl doped ice at -35
o
C,
the ion concentration in the ice is approximately 2 μM near saturation as shown on Figure
3.6.3.5. Figure 3.6.3.6 summarizes the variations of ion saturation concentrations in the ice from
different ionic species as reference to pure ice. As expected, the ion concentration decrease with
decreasing temperature which is similar to the behavior of defects. Ice doped with KCl, NaCl,
and KI exhibited similar concentrations near 2 μM. The corresponding ionic diffusion
coefficient as estimated from the Nernst-Einstein equation is about 10
-5
cm
2
/sec which is within
expectations of realistic values. Ice doped with MgSO
4
exhibited the lowest concentration of
approximately 0.04 μM. Again, the highest ion concentration is ice doped with KOH and CaCl
2

at nearly 2 orders of magnitude higher, but is most likely due to liquid-like conditions with free
ion motion.
By asserting the classical interpretation that Bjerrum protonic defects dominate the
observed bulk ice conductivity, the ion saturation concentration is equivalent to the effective
defect concentration. Given the effective partial charge of Bjerrum defects is approximately
0.38e (Petrenko and Whitworth 1999, pg. 154), this corresponds to approximately 1.36 x 10
15

defects/gm of ice which agrees with previous findings.
48,49
Therefore, the corresponding
mobility in KCl-doped ice is approximately 4.0 x 10
-7
m
2
/V-sec and is also in good agreement
with other studies.
50
This defect mobility is substantially lower than that of free ion motion
observed in liquids, which appear to substantiate Bjerrum type defects as a dominate factor in
contributing to the observed bulk conductivity of ice. However it is important to reiterate that
this defect concentration in ice is really facilitated by ion absorption. In addition, a more
accurate understanding of ionic conduction in the ice lattice need to be further addressed in order
to develop realistic mechanisms for defect motion in the ice. Chapter 4 attempts to ascertain
such details on the microscopic dynamics of ice though atomistic simulations.

101

Figure 3.6.3.5: Ion concentration plot for KCl doped ice at T = -35
o
C. The conductivities of
100 mM and 1 M KCl doped ice are at the ion saturation concentration of about 2 μM. This ice
ion concentration is estimated using the calibration measurements of conductivity on liquid
electrolytes.

Figure 3.6.3.6: Summary of ion saturation concentration in ice for various dopant types as
reference to pure ice.

102

3.7 Summary & Conclusions
The investigation of pure and ion doped ice by EIS is a viable and robust technique for
in-situ characterization on the unique electrical properties of ice as emanated from its molecular
structure. To demonstrate this, we employed a novel point electrode configuration and a housing
apparatus capable of growing of single-crystalline ice and conducting in-situ measurements.
Measurement results on pure ice are in excellent agreement with measurements from baseline
parallel plate electrodes. From equivalent circuit modeling, the Debye relaxation of our point
electrode configuration was determined to be 0.579 eV as oppose to 0.529 eV from the plate
electrodes. Resonance is also detected with the point electrode configuration at the onset of the
Debye relaxation which is unique to our measurement results. The detected resonance is caused
by the radial orientation of the applied electric field and is responsible for the measurement
differences as compared to the parallel plate electrodes.
In order to investigate the effect of ionic species on electrical properties of ice,
measurements of ice doped with KCl, NaCl, KI, KOH, and CaCl
2
were performed over a large
temperature range down to -65
o
C. Ice doped with KOH and CaCl
2
, however, appear to inhibit
the freezing process down to as low as 65
o
C resulting in liquid-like behavior or amorphous-like
ice structure. In all these cases, conductivity exceeded that of pure liquid water at values of
about 10
-4
S/cm indicating free ion motion. On other cases, ice from KCl, NaCl, and KI
exhibited very distinct ice signatures of a single Debye relaxation and resonance along with the
equilibrium polarization states as described by ε
s
and ε
∞
in the relative permittivity spectra. The
main significant difference is the Debye relaxation which occurs at significantly higher
frequencies than that of pure ice. The measurement signatures of ice from MgSO
4
are the most
similar to that of pure ice suggesting the lowest ion concentration suggesting almost all ions
ejected from the ice structure during freezing. This was further substantiated by visualization of
pure and doped ice grown in polystyrene cup, where it was observed that ice nucleation occurs at
the outer wall towards the center. The process of ion-pairing appears to exist along the same
nucleation pathway causing aggregation of cracks and fractures to form at the center of the EPC.
Thus, it is concluded that ionic type play an important role in the freezing mechanism of ice in
which they can migrate freely, trapped in ice cavities, or totally ejected causing salt
recombination and phase separation observed by EIS measurements.

103

The time constants for relaxation and resonance were determined using an equivalent
circuit which accurately predicts both relaxation and resonance behavior. As expected, the
Debye relaxation and resonance time constants exhibited very similar trends. Typically, ice with
halogenated salts exhibited the lowest time constants while ice from MgSO
4
has the highest
values. From classical theories, all solvated ions in the ice lattice are typically trapped in the ice
cavities and are immobile. However, these absorbed ions also facilitate defect formation and
migration which ultimately resulted in lower relaxation time constants. On the other hand, ions
from KOH and CaCl
2
exhibited higher activation energies than that of halogenated salts but
conductivity is substantially higher by more than 2 orders of magnitude. For such ionic species,
the freezing process and ice crystallization appear inhibited resulting to amorphous type ice
containing substantially more free-ion motion as the main factor in carrying current. Finally, in
the case MgSO
4
the time constants are the closest to that of pure ice with the highest
corresponding activation energy which also indicate both low ice concentration and
corresponding negligible changes in defect formation. Therefore, ion type is presumed to play a
critical role during the freezing process as indicated by the Debye relaxation and resonance times
and their corresponding activation energies. Though the activation energies for relaxation and
resonance are vital factors in facilitating defect formation and migration, the most dominating
factor remains ion and/or defect concentration.
From conductivity measurements on doped ice, quantification of Bjerrum defect
concentration and mobility in the non-polarized frequency region can be estimated up to the ion
saturation concentration. The mechanism for conductivity appears to vary greatly due to ionic
species as well. Despite highest measured conductivity, ice doped with KOH also exhibited by
highest activation energy of 0.495 eV, and suggests free-ion motion plays a dominate role as
oppose to solvated halogenated salts such as KI, NaCl, and KCl which exhibited substantially
lower activation energies despites despite their anticipated lower ion motion. In the case of KCl
doped ice, the estimated defect concentration at the ion saturation concentration was also
estimated at 2 μM with a corresponding ionic diffusion coefficient of about 10
-5
cm
2
/sec. An
equivalent defect concentration is approximately 1.36 x 10
15
defects/gm of ice. This corresponds
to a mobility of about 4.0 x 10
-7
m
2
/V-sec which is in excellent agreement with previous studies.

104

Chapter 4: Molecular Dynamic Simulations of Ice

4.1 Literature Review of Ice Simulations
A substantial body of research exists on ice dating back to the late 1920’s. From over 80
years of research, we witnessed an evolution of ice research from theoretical formulations to
experimentations to MD modeling and simulation as illustrated in Figure 4.1.1. However, the
most recent atomistic simulations of ice focus on either ice with ion or ice with an applied
electric field but not both. Our current atomistic study of ice will further advance our
understanding of the molecular physics of ice by incorporating both ion and an applied electric
field in the simulation system.

Figure 4.1.1: Summary on historical timeline of ice research.

Atomistic simulations of water/ice have expanded tremendously due to the recent
advances in computer power. In regards to simulations of water/ice, research over the past
couple of decades, molecular dynamic (MD) and monte carlo (MC) simulation have increased in
systems size along with capabilities to include simple ions and even large molecules. This
enabled the investigation of many types of studies, including determination of phase diagrams,
freezing/melting mechanisms, transport processes, and even direct calculation of electrical
properties desirable in this investigation. The following literature review literature review
1920 1930 1940 1950 1960 1970 1980 1990 2000 2010 2020
Ice Research Timeline
Theoretical Ice Formulations:
Bernal & Fowler, 1933.
Pauling, 1935.
Bjerrum, 1952.
Ice Experimental Research:
Auty & Cole, 1952.
Camplin, Glen, and Paren, 1978.
Kawada, 1987 & 1997.
Itoh, Kawawmura, Hondoh, Mae, 1996
MD Bulk Ice:
Buch, Sandler, Sadlej, 1998.
Smith, Haymet, 2004.
Vega, McBride, Sanz , Abascal, 2005
MD Ice + Ion:
Ikeda-Fukazawa, Horikawa, Hondoh, 2002.
Smith, Byrk, Haymet, 2005.
MD Ice + E-field:
Vegiri , Schevkunov, 2001.
Wei, Zhong, Su-Yi, 2005.

105

focuses on critical investigations of ice and analysis techniques used on results from atomistic
simulations.
Analysis of atomistic simulations of ice systems followed efforts of Bernal and Fowler
(1933), Pauling (1935), and Bjerrum (1952), which formulated the theoretical framework of ice
from a molecular perspective. Auty and Cole (1952) provided one of the most complete
experimental investigations on the electrical properties of ice. The experimental investigation
focusing specifically on ice transport mechanisms was conducted by Goto et al. on measuring the
diffusion coefficients of self-interstitials.
51
The only experimental-based measurement of the
radial distribution function (RDF) in water and ice using neutron diffraction was provided by
Soper.
52
Together, these seminal investigations formed the basis of approach for analysis on the
structure and transport properties from atomistic simulation such as MD/MC in ice and/or water
systems appropriate for interpreting macroscopic electrical properties.
To enable atomistic studies, important initial works on ice emphasize more toward the
development of a valid ice lattice by incorporating both theoretical and experimental insights.
Abascal et al. compared the physical properties determined from MD simulation on all the water
interaction potentials, and found results from TIP4P water model is in best agreement with
experiment.
53
This study also led to the development of TIP4P/ice and TIP4P/2005 water
models specially designed for simulations of ice. A thorough investigation of structural stability
and densities of SPC/E, TIP4P, and TIP5P potential models on ice was performed by Vega et al.
comparing the radial distribution functions of O-O from NPT ensemble simulations.
54
It was
concluded that proton disorder contributes significantly to structure differences from potential
model, and the formation of doublet H atoms (D-defect) is feasible from order-disorder
configurations. To establish the potentially favorable positions for ion solubility, Smith and
Haymet performed an energy contour profile of Na
+
and Cl
-
atoms in ice Ih lattice.
55
They
estimated the minimum potential energy of -9.0 kcal/mole corresponds to the hexametric cavities
of the lattice. The most feasible method of determining energetic favorability of ion locations in
the ice lattice was described by the work of Feibelman.
56
From First-principles calculations, this
study found Na
+
and Cl
-
ion hydration in ice is more energetically favorable by substitution of
water molecules in the lattice as oppose to interstitial ion solvation. The conclusions from this

106

study validate a more simplistic technique of ion insertion for water molecules in simulation of
large-scale ice/ion lattice without conducting a painstaking energy minimization analysis.
The majority of early atomistic studies on ice centers on the structure and dynamics of
orientational defects. These MD studies on ice essentially focused on contributions of proton
disorder to support observation from critical experimental studies.
57
For example, Itoh et al.
investigated the librational, vibration, and translation degrees of freedom on atomic motion on
ice as a basis for interpreting experimental neutron scattering and power spectra.
58
Podeszwa
and Buch were one of the first investigators to study the role of such orientational defects on
phenomenological processes such as dielectric relaxation and electrical conductivity (Podeszwa
et al., 1999). Their study included the manual insertion of an actual Bjerrum defect pair of D
(excess OH) and L (mission OH) defect by rotation of the water molecule prior to performing the
MD simulation. The analysis included the examination of the trajectories of the simulation and
computation of the lattice energy with and without the presence of defect to assert the conclusion
of defect migration caused by vibrational coincidences. The seminal studies by Ikeda-Fukazawa
on the diffusional mechanism of ice Ih took a step further by including the presence of ions on
defect formation and migration.
59,60
To incorporate the ions, they used the method of calculating
the energy contour profile near the active sites to determine the energy minima. Their
comprehensive analysis included the trajectory plots for visualization of ionic motion and
square-mean-displacement (MSD) to calculate the diffusion coefficient at various temperatures.
Diffusion coefficients were found to be in the range of 10
-9
to 10
-11
m
2
/sec depending on the
molecular size of the migrating ions or molecules. The activation energy of the diffusion
coefficient was determined using the Arrhenius equation. Their results concluded that there exist
two unique mechanisms of migration: 1) vibrational motion of interstitial molecules for smaller
molecules, and 2) the breaking of hydrogen bonds in the ice lattice suggesting Grotthuss hopping
mechanism at work for larger atoms or molecules. Their conclusion provided a basic
understanding of the ionic diffusion mechanism in the lattice.
Ionic hydration in the ice Ih lattice represents another substantial body of work. One of
the first comprehensive work on this was performed by Kalko et al.
61
This work focused on the
dielectric contributions of divalent and trivalent cation and anions to the free energy using simple
ions of Na
+
, Ca
+2
, and Cl
-
by using the parameter of the covalent radius of conducting hard

107

sphere as predicted by the Born model. The study calculated the coordination number from RDF
RDF. To analyze effect of ions on water molecules, the study calculated residence time by
determining the number of molecules that lie initially within the 1
st
coordination shell and are
still there after an elapsed amount of time. The effect of system size and potential model on
solvation free energies of methane was studies by Lyubartsev et al. using AMBER MD
package.
62
The significance of this work mainly found that the solvation of organic materials is
impacted by the presence of ion hydration such as NaCl. The work of Byrk and Haymet
employed a larger system of 2304 SPC/E water molecules to study the Na
+
and Cl
-
solvation
mechanism in surface ice.
63
This work analyzed the single-ion density profile to show ion
density at different simulation times to suggest oscillation of ions occur prior to diffusion. A
running coordination number was also calculated to probe the molecular environment around the
ion. Results show the Cl
-
ion diffuse more readily than the Na
+
ion, mainly attributable to ionic
size and concentration. A similar study by the same research group also investigated the Na
+
and
Cl
-
solvation across the water/ice interface.
64
This work included a thorough structural analysis
by RDF of Cl-O and C-H atoms at temperatures of 298
o
K (water) and 150
o
K (ice). RDF analysis
indicated melting conditions near the ions while the away from the ions the ice structure is
stable. The coordination number also revealed that the water and ice have unique structure near
ions. Other MD studies on the melting behavior at the water/ice interface also revealed similar
conclusions attributable to ionic diffusion across the water/ice interface.
65,66

Non-equilibrium molecular dynamic (NEMD) studies involving the addition of an
applied electric field have enabled numerous new types of investigations applicable to ice. In
liquid water, the application of an electric field induces various degrees of orientation and
translational order causing ice-like arrangements and a process called electro-freezing.
67
The
applied electric field strength of such MD studies typically exceeded 7 V/nm to achieve the
desired effect. Evan et al. observed this phenomenon with applied electric field strength of 7
V/nm in the investigation of thermal transport characteristics.
68
The study used a combination of
RDF, structural factor, mean-dipole moments to monitor a crystalline arrangement like ice.
Another interesting study by Yan and Patey also applied an electric field gradient in a MD study
to investigate ice nucleation sites relative to the electric field penetration depth.
69
For optimal
accuracy, the water model TIP4P/ice was chosen with the dipole moment of 2.43D. Their
analysis also included RDF and MSD at various temperatures to achieve a conclusion that the

108

nucleation process is independent electric field strength. Structural analysis in the presence
electric field is typically performed by examining the H-bonding network. By using a criteria of
a 6-bond hydrogen ring for liquid conditions Hyun Jung et al. observed in TIP4P water that the
applied electric field in water stabilized liquid phase creating a crystal-like behavior.
70
The
diffusion coefficient was found to significantly decrease with the applied field. In a similar
structural study, Wei et al.
71
employed analysis of RDF O-O and O-H atoms in addition a
quantity called a hydrogen bond correlation function to account for the time-dependent bond
variable which equal to 1 if pair of water is hydrogen bonded at time t and 0 otherwise. This
method describes the structural relaxation of hydrogen bonds and the associated relaxation time
can be interpreted as time scale of reorganization of hydrogen bonds through the calculation of
the fraction of water molecules engage in hydrogen-bonds (range from 1-6), and the average
number of H-bond per water molecule. This study also reported anisotropy effects from
calculating diffusion coefficients along the X,Y,Z-axis, and also concluded that the diffusion
coefficient decreased by a factor of 11 times compared to that of no field. Another seminal
investigation using NEMD was performed by Kadota et al.
72
This focused on the dehydration
process of supersaturated NaCl solution under the influence of an electric field with analysis
focus on the ionic motions at the solid-liquid interface. The density profiles and the dynamic
diffusion coefficient of water molecules as a function of distance from the crystal interface were
estimated by characterizing such interfaces. Dehydration energy for crystal growth is also
determined, with 410 kJ/mole for Na
+
ion and 317 kJ/mole for Cl
-
ion. It was found that the
structure and dynamics of the interfaces were different from those of the bulk solution as a result
of both the electric fields and the surface charge; as the electric field strength increased, the Na+
ions broke out of the hydration structure but contradictory to ice simulations the overall
hydration structure in the bulk solution is not affected. Overall, such NEMD studies can provide
further insights from unique perspectives on the dynamic structure for crystallization and ionic
solvation applicable for both water ice phases.
Perhaps the most thorough MD studies on the structural properties with an applied
electric field using NEMD was performed by Vegiri and Schevkunov.
73,74
Their studies included
an application of various electric field strength in the order of 10
7
– 10
8
V/cm at T = 200K using
TIP4P water model. Standard energy profiles and trajectory plots were also employed in this
study. However, the in-depth analysis also included the normalized RDF of O-O atoms along

109

with the corresponding oxygen atom coordination number, total dipole moment fluctuations,
MSD, and translational diffusion. They reported that the highest applied field strength induces
reorientation leading to phase change with square stacked ice molecules arranged
perpendicularly to the field direction associated with electrofreezing. However, it was also
found that weak fields cause a transition from the solid to the liquid state. Such studies in
general supports an overall conclusion that the effect of applied field strength on structural and
dynamics changes is complex and non-monotonic ways which have associated impact on the
chemical reactivity and the solvation capacity of ice.
Few direct atomistic simulation studies exits on the electrical properties of ice such as
dielectric constant and ionic conductivity. Rick and Haymet conducted the only investigations
on the static dielectric constant, ε
s
, due to proton disorder using MC simulations.
75
The values of
ε
s
were determined by the equation over a range of temperatures from 10 to 273
o
K. The optical
dielectric constant, ε
∞
, used in this study was 1.592. Results from this study strongly supported
actual experimental measurements of ε
s
performed by Auty and Cole (1952). Though not on ice
directly, the computation of ionic conductivity from MD simulations was performed on other
material systems of solid state which can serve as a basis for ice electrical properties analysis.
Tang et al. investigated oxygen mobility in yttria stabilized zirconia at temperature in excess of
1000
o
K with an applied electric field in the order of 10
7
to 10
8
V/m.
76
Noritake at al. compared
the conductivity calculations from MD simulations in order to investigate conductivity
measurements in Na
2
O
.
3SiO
2
melts under high pressure.
77
This study determined average
diffusion coefficients from mean-square-deviation (MSD) data. From the diffusion coefficient,
the conductivity was calculated by using a form of the Nernst-Einstein relation. The activation
energy was also determined using the Arrhenius equation. A more general direct method of
calculating the electrical conductivity contributions of cations and anions was employed by
Koishi and Fujikawa.
78
This MD study also added an electric field corresponding to an applied
force of 0.02 to 1.60 x 10
-9
N to determine current yielded by the external force. The electrical
conductivity of the cations and anions which was referred to as the partial electrical conductivity
is then estimated using the known volume of the simulation cell. Therefore, the total electrical
conductivity is the sum of each of the partial conductivities (σ=σ
+
+ σ
-
) from both anions and
cations.

110

4.2 Theory and Simulation Details
This chapter explores the microscopic dynamics of ice systems through “computer
experiments” as pioneered more than 40 years ago.
79,80,81
From the two general categories of
atomistic simulations, namely molecular dynamics (MD) and Monte Carlo (MC), it is the MD
technique that has the advantages of providing dynamic information of the system though
solving the equations of motion directly and follow the evolution of the positions and velocities
of the particles. These equations of motion that allow the atoms/particles to move as a function
of time are governed by Newton’s second law
. ,..., 1 ,
2
2
N i F
t
r
m
i
i
i
 


4.2.1
For an N-particle system there are hence three differential equations for each particle i of mass
m
i
, located at point r
i
= (x
i
, y
i
, z
i
) in Cartesian coordinate space. Given the known forces, the
motion of the system as a function of time, r
i
(t), can be achieved by solving the 3N second-order
differential equation of equation 4.2.1 coupled with 6N initial conditions of the positions, r
i
(t=0),
and the velocities, v
i
(t=0).
From the framework of classical mechanics where the forces depend only on the
coordinates and not explicitly on time, the forces can be derived as the gradient of the potential
function
) ,..., (
1 N r i
r r U F
i
  4.2.2
where the subscript emphasized that the gradient is to be taken with respect to the position of
particle, i. Hence, the forces are conserved, since the work done by moving a particle between
two points is independent of the path taken. For simulations with ion doped ice, the main
contributions to the pair potential forces are the non-bond interactions from the strong repulsive
force at very short distances due to electron cloud overlap (the Pauli principle) and the attractive
force at moderate distances due to correlation between the fluctuating multipoles (London
dispersion) as famously described by the Lennard-Jones potential, coupled with the electrostatic
forces from ions (Coulomb charge-charge interactions). The total pair potential is described by

111

ij
j i
ij ij
ij
r
q q
r r
r E
) (
4 ) (
6 12




























 
 4.2.3
where r
ij
is the inter-particle separation radius, ε and σ are standard Lennard-Jones (LJ)
parameters, and q
i
and q
j
are the charges of two interacting particles i and j, respectively. For
interaction between unlike particles, LJ parameters are calculated from self-interaction
parameters employing Lorentz-Berthelot combining rules:
j i ij
j i ij
  
  

  ) (
2
1
4.2.4
We assume bonds are rigid and the total potential is pairwise additive.
One of the key technical challenges of this investigation is the selection of an appropriate
transferable intermolecular potential (TIP) force-field that would generate stable, correct ice
structure. Typical TIP water models appropriate for ice simulations include three or four-points
models as summarized in Table 4.2.1. For our study, we employed only the four-point water
models (TIP4P) which contain a bisector point M as shown in Figure 4.2.1. We developed ice
lattices comprising of several types of TIP4P molecules using a unit cell in a cubic arrangements.
The unit cell types containing 8 water molecules are shown in Figure 4.2.2. The O-O separation
distance is 2.75 Å and the H-O-H dihedral angle is 109
o
. The ionic radius for K
+
ion is 2.72 Å,
and for the large Cl
-
ion is 3.79 Å. From a periodic expansion technique, supercells of various
sizes (1600 and 6144 water molecules) were developed from the unit cell for production
simulation runs.

112

Table 4.2.1: Physical Parameters for Ions and Water 3-Point and 4-Point Model.

The blue highlighted force-field parameters were employed in this study.

Figure 4.2.1: Illustration of four point water models for atomistic simulation. The M point is not
depicted accurately due to very short OM radius.

Parameter TIPS SPC TIP3P SPC/E BF TIPS2 TIP4P
TIP4P-
Ew
TIP4P/
Ice
TIP4P/
2005
K
Atom
Cl
Atom
r(OH), Å 0.957 1 0.96 1 0.96 0.96 0.96 0.96 0.96 0.96 --- ---
HOH, deg 104.5 109.5 104.5 109.5 106.0 104.5 104.5 104.5 104.5 104.5 --- ---
r(OM), Å --- --- --- --- 0.15 0.15 0.15 0.13 0.16 0.16 --- ---
q(O or M) −0.80 -0.82 -0.83 −0.85 −0.98 −1.07 −1.04 −1.05 −1.18 −1.11 --- ---
q(H) 0.40 0.41 0.42 0.42 0.49 0.54 0.52 0.52 0.59 0.56 --- ---
r(Ion), Å --- --- --- --- --- --- --- --- --- --- 2.72 3.79
q(Ion) --- --- --- --- --- --- --- --- --- --- +1.0 -1.0
Molar Mass,
gm/mole
18.0 18.0 18.0 18.0 18.0 18.0 18.0 18.0 18.0 18.0 39.1 35.4

113

Figure 4.2.2: Ice unit cells of a cubic arrangement.

In order to create ice lattices comprising of four point water models (e.g., TIP4P, TIP4P-
Ew, TIP4P/Ice), an algorithm was developed to compute the xyz-coordinates of the point (M)
along the water molecule bisector. The algorithm involves solving a set of equations from vector
geometry. Equation 4.2.5 is used to find the equation of the plane corresponding to each water
molecule of the form
D cM bM aM
z y x
   4.2.5
where the coefficients a, b, c, and D are solved by definition of cross product and known vector
coordinates of the two OH vectors in the plane of the water molecules. Equation 4.2.6 is the use
of the inner dot product equation for angle H-O-M such


Cos OM OH O M O H O M O H O M O H
OM OH
OM OH
Cos
z z z z y y y y x x x x
        


) )( ( ) )( ( ) )( (
4.2.6
where OH and OM are vectors, |OH| and |OM| are magnitudes of the vectors OH and OM ,
and θ is the bisector angle (e.g., 15
o
for TIP4P model) of angle HOH. Equation 4.2.7 is the
definition of the perpendicular dot product between vectors
2 1
H H and OM such that

114

0 ) )( ( ) )( ( ) )( (
0
1 2 1 2 1 2
2 1
        
 
z z z z y y y y x x x x
O M H H O M H H O M H H
OM H H
4.2.7

Thus, solving equations 4.2.5, 4.2.6, and 4.2.7 by matrix Gaussian elimination techniques yields
the coordinates of M
x
, M
y
, and M
z
for each water molecule in the entire lattice network. All atom
coordinates of the lattice are outputted to a single file in pdb (protein database) format. The data
of the pdb file was subsequently converted to corresponding formatted inputs files of initial
coordinate compatible for MD simulation using AMBER or GROMACS. For comparison
purposes, ice lattices of all the four point water models (TIP4P, TIP4P-Ew, TIP4P/Ice, and
TIP4P/2005) were created for MD simulation. Figure 4.2.3 illustrates the ice lattice developed
from different water model containing 1600 water molecules. All four ice lattices (A-D) are
identical with exceptions of the individual water molecules making up the entire ice lattice. The
ice lattice satisfies to the ice rules of hexagonal Ih ice. Figure 4.2.4 shows the ice lattice
containing 6144 water molecules using TIP4P model for simulations containing ions and the
influence of an electric field.

(A) (B)

115

Figure 4.2.3: Ice lattice using different water models: (A) TIP4P model, (B) TIP4P/Ice, (C)
TIP4P/2005 Model, (D) TIP4P-Ew. All ice lattices contain 1600 water molecules with box
dimensions of 39.3 Å x 38.9 Å x 36.7 Å.

Figure 4.2.4: Example of ice lattice containing 6144 water molecules with box dimensions of
53.9 Å x 62.2 Å x 58.7 Å.

To develop ice/ion lattices for simulation, we employ a substitutional technique of
randomly replacing water molecules for ions as supported by the work of Feibelman (2007).
This is especially appropriate for the larger Cl
-
ions relative to the space available within the ice
cavities. Three unique initial ion/ice lattices of the same concentration were developed in order
to study the effects of initial condition of ion locations in order to ensure repeatability and
provided statistics from the results. The possibility of creating a false migration pathway in the
(C)
(D)

116

ice lattice though creation of vacancy defects is minimized by applying a minimum of 2.5 Å
separation radius between ions. In the 3 initial ion configurations, configuration #2 has radius
separation greater than 10 Å between ions. The lattice dimensions of these lattices are 53.9 Å x
62.2 Å x 58.7 Å. We chose 30 ion pairs (60 ions total) corresponding to a concentration of
approximately 0.253 M in order to simulate conditions of briny concentrations similar to the
experimental systems. Figure 4.2.5 illustrates the 3 initial configurations in a 3D ice lattice
containing ions.

Figure 4.2.5: Initial configuration of a 3D ice lattice containing 30 ion pairs of K
+
and Cl
-
ions
(60 ions total). The large while dots are K
+
ion and the large blue dots are Cl
-
ions. (Left) initial
ion configuration #1, (Middle) initial ion configuration #2, (Right) initial ion configuration #3.
Initial ion configuration #2 is minimum ion separation radius of 10 Å.
Both open source MD packages of Amber (version 10.0) and Gromacs (version 4.5.5)
were used to perform MD simulations. In general, Amber simulations for performed for
purposes of initial MD setup and validation of available TIP4P water force-field models. Amber
was also used for conduct smaller simulation of less than 2000 molecules. Upon successful
validation of such smaller simulations, Gromacs was employed to perform large productions
simulations containing more than 6144 molecules depending mainly due to its capabilities of
adding an external electric field. The complete time-dependent electric field pulse equation
82
in
Gromacs is of the form
  ) ( cos
2
) (
exp ) (
2
2
o
o
o
t t
t t
E t E 





 
  

4.2.8

117

where ω is the angular frequency, σ is the width of the pulse, t
o
is time at pulse maximum, and
E
o
is the amplitude of the applied electric field. Gromacs also has an option to convert the full
pulse equation to the continuous sinusoidal wave function of the form
) cos( ) ( t E t E
o
  4.2.9
and is deemed suitable to represent the applied electric field from a impedance spectroscopy
experiment. Using equation 4.2.9, we applied electric field strength amplitude (E
o
) of 1 V/cm
with angular frequency of 1 MHz which corresponds to a nearly constant electric field over a 1
μ-sec simulations duration. These electric field parameters for simulation are virtually identical
to the experimental conditions of the applied electric field and, therefore, offer direct comparison
of simulation results to experimental results.
Both Amber and Gromacs MD setup are very similar. All simulations apply SHAKE
algorithm for rigid water molecules to constrain hydrogen covalent bonds, and the integration of
the dynamic equations was done using the leapfrog algorithm. Both the short-range electrostatic
interaction and the Van der Waals interaction are cut-off at 12 Å. The main difference between
the Amber and Gromacs simulations is the thermostat. Amber thermostat employs the Anderson
temperature coupling scheme in which distribution of randomize velocities correspond to
assigned temperature. In Gromacs simulations, both Berendsen and Nose-Hoover thermostats
were used. The barostat is an isotropic position scaling algorithm. The simulations were
performed with a time step of less than 0.002 ps depending on the simulation size. The
simulation temperatures of the production runs incorporating ions include both 200K and 220K.
In order to thermalize the system, a temperature step of 5K was employed to heat or to cool the
system for 10 ns for each temperature production run. An additional 20 ns initial simulation was
conducted for purposes of energy minimization and to establish thermal equilibrium using NVE
ensemble without electric field. System relaxation was also performed after ion insertion for 40
ns. Production MD runs was conducted using NpT ensemble of 20 ns simulation duration with
the addition of an electric field. On all production runs, each MD restart run employs the
previous coordinate configuration as its initial coordination configuration until a simulation
duration of at least 1 μ-sec was achieved. Each production run took approximately 40 days of
actual CPU time using 96 CPU nodes of the HPCC cluster at USC.

118

4.3 Results and Discussion
As illustrated in by various previous studies, there are numerous statistical formulations
available to organize and to analyze the simulation data in an effort to probe the physical,
structural, and transport dynamics of ice. Our study on ice mainly emphasizes the role of ions as
related to the observed electrical properties of ice through calculation of radial distribution
functions (RDF), mean-square-distribution (MSD), and diffusion coefficients. To evaluate the
lattice stability and quality, we characterized the system energetics and calculated some critical
physical properties for comparison with known values. The RDF was also used to ascertain
structural stability in terms of oxygen distribution and to provide critical information on the
structural dynamics of the ice lattice as impacted by both ions and the externally applied electric
field. Details on ionic migration mechanisms were characterized using individual ionic
displacement profiles. The MSD and the corresponding diffusion coefficients enable the
computation of the bulk electrical property of ionic conductivity for direct comparison with
experimental results. To examine the local structure surrounding the ions and to visualize the
ionic motion within the ice lattice, we also employed various forms of trajectory plots.

4.3.1 Kinetic Energy, Potential Energy, and Density
Prior to performing large ice lattice simulations, MD simulation of smaller ice lattices
comprising of different TIP4P water models were investigated for optimal structural stability.
These studies runs were conducted in AMBER and comprise of only 1600 water molecules. As
illustrated in Figure 4.3.1.1, after an initial rapid energy rise due to thermal equilibrium all four
types of ice lattices achieved equilibrium after 20 ps simulation time and exhibited excellent
equilibrium kinetic and potential energies at the end of the 100 ps simulation time. Unlike the
kinetic energy, the potential energy is dependent on the water model type. Table 4.3.1.1
summarizes the molar energies of ice obtained from various water models. The highest
equilibrium total energies were exhibited by lattice containing the TIP4P and TIP4P-Ew water
models at -46.4 kJ/mole and -45.9 kJ/mole, respectively, and compares well with the reported
experimental value of -47.3 kJ/mole (Petrenko and Whitworth 1999, pg. 29). The ice lattice with
TIP4P/ice model exhibited approximately 25% lower equilibrium potential energy as compared
to the TIP4P model at -53.5 kJ/mole.

119

Large ice lattice simulations using the TIP4P model were also employed to minimize the
effects of periodic boundary conditions (PBC). The dependence of potential and kinetic energy
on temperature for the pure ice system containing 6144 water molecules are shown in Figure
4.3.1.2. Both the kinetic and potential energies of the system increase with temperature as
expected, with equilibrium energies nearly identical to simulations using 1600 water molecules.
The kinetic energy at 200K is 5.0 kJ/mole and increased about 0.60 kJ/mole at 220K. The
observed potential energy increased from 200K to 220K is approximately 0.66 kJ/mole at -46.4
KJ/mole and -45.3 KJ/mole, respectively. Due to negligible pressure and volume changes of
condensed materials, the heat capacity can be determined from the total internal energy
difference from the 200K and 220K simulation temperatures. The simulation heat capacity is
approximately 57.9 J/K-mole which is in good agreement with actual reported values for ice
(Petrenko and Whitworth 1999, pg. 41).

120

Figure 4.3.1.1: Total energy profile comparison of all four-point water models containing N=
1600 water molecules: (A) Potential energy, (B) Kinetic energy. Simulation performed in
AMBER. Simulation duration = 100 ps.

(A)
(B)

121

Figure 4.3.1.2: Energy profile comparison of N= 6144 water molecule system of pure ice using
TIP4P/Ice water model at different temperatures: (A) Potential energy, (B) Kinetic energy.
Simulation performed in Gromacs. Simulation duration = 40 ps.

With the addition of electric field, both the potential and kinetic energy exhibited no
appreciable difference from 200K to 220K as compared to the simulations without the presence
(A)
(B)

122

of the electric field. As shown in Figure 4.3.1.3a, the kinetic energy increased with temperature,
with the average kinetic energy of 5.0 +/- 0.03 kJ/mole at 200K and 5.5 +/- 0.03 kJ/mole at
220K. As expected, the potential energy also increased by the same amount as compared to the
simulation without the electric field as shown in Figure 4.3.1.3b. The standard deviation also
remains low at +/- 0.03 kJ/mole. Therefore, the change in total system energy is negligible in the
presence of the electric field, but an increase of about 1.0 KJ/mole is observed from 200K to
220K.
Unlike the electric field, the effect of ions on system energies is more significant. Over
the 1 μ-sec simulations time, both the potential and kinetic energies are relatively stable as
shown in Figure 4.3.1.4a,b. Like pure ice, the potential energy increased an average of 0.50
kJ/mole with increase in temperature from 200K to 220K. Under the influence of the electric
field, the estimated heat capacity in these simulations systems is approximately 51.4 J/K-mole.
In addition, the ice lattices containing K
+
and Cl
-
ions exhibited higher standard deviations in
potential energy, by as high as a factor of 2 as compared to that of pure ice. As tabulated in
Table 4.3.1.2, the standard deviation of both simulation systems containing ions is about +/- 0.06
but is independent of both temperatures and the presence of the electric field. Therefore, these
observed potential energy differences suggest the presence of the ions in the ice lattice that have
a more dramatic effect over the influence of the electric field.

123

Figure 4.3.1.3: Energy profiles comparison of pure ice system containing 6144 water
molecule under the influence of the electric field: (A) Potential energy, (B) Kinetic
energy.

(A)
(B)

124

Figure 4.3.1.4: Total energy profile comparison of ice + ion system containing 6084
water molecules and 60 ions under the influence of the electric field: (A) Potential
energy, (B) Kinetic energy.

(B)
(A)

125

By similar comparisons, all TIP4P water models resulted in a lattice density consistent
with real ice. The lowest density was exhibited by TIP4P/ice model shown on Figure 4.3.1.5
with an average density of 921.1 Kg/m
3
which is in good agreement with the reported
experimental value of ice density of 920.0 Kg/m
3
at 260K (Petrenko and Whitworth 1999, pg.
29). The TIP4P-Ew, TIP4P/2005, and TIP4P water models exhibited average densities of 923.1
Kg/m
3
, 931.8 Kg/m
3
, 946.3 Kg/m
3
, respectively. The simulated densities at 200K and 220K
from larger lattices are shown in Figure 4.3.1.6a,b. The density average and the standard
deviation comparing the effects of temperature, ions, and the influence of the electric field are
also summarized in Table 4.3.1.2. The density at 200K is 945.7 +/- 0.09 Kg/m
3
and decreased to
942.0 +/- 0.02 Kg/m
3
at 220K for the pure ice system under the influence of the electric field.
The differences appear insignificant as the average density difference is less than 0.5%
comparing 200K and 220K simulation results. In the ice and ion system under the influence of
the electric field, the average density increased to approximately 959.3 Kg/m
3
as compare to the
pure ice system. Similar to the pure ice system, however, density along the standard deviation
also decreased with increase in temperature from 200K to 220K. By no means definitive, this
suggests the addition of ions in ice may cause structural or packing changes. Overall, the
simulated ice density is well within the criteria of realistic values even with the added ions.

Figure 4.3.1.5: Density profile comparison of TIP4P, TIP4P/ice, TIP4P-Ew, and
TIP4P/2005. The ice system comprise of 1600 water molecules with no applied electric
field.

126

Figure 4.3.1.6a: Density profile comparison of 6144 molecule pure ice system under the
influence of the electric field.

Figure 4.3.1.6b: Density energy profile comparison of ice + ion system under the
influence of the electric field.

127

Table 4.3.1.1: Comparison of average total potential and kinetic energies and average density
from Various Water Models at T = 200K. All simulation systems contain 1600 water molecules.

Table 4.3.1.2: Summary of System Molar, Potential, and Kinetic Energies and Densities of
various simulation systems at T = 200K and T = 220K. Pure ice simulation systems contain
N=6144 molecules. Ice + ion systems contain N=6084 water molecules and 60 ions.

Water Model
Potential Energy
(kJ/Mole)
Kinetic Energy
(kJ/Mole)
Avg Total Molar
Energy (kJ/Mole)
Avg Density
(Kg/m
3
)
TIP4P -51.4 5.0 -46.4 946.3
TIP4Pew -50.8 5.0 -45.9 923.1
TIP4P-Ice -65.3 5.0 -60.3 921.1
TIP4P/2005 -58.5 5.0 -53.5 931.8
Temperature
Pot Energy
(kJ/Mole)
Kinetic Energy
(kJ/Mole)
Total Energy
(kJ/Mole)
Density
(Kg/m
3
)
200 -51.4 +/- 0.03 4.9 +/- 0.03 -46.4 +/- 0.01 945.7 +/- 0.09
220 -50.7 +/- 0.03 5.5 +/- 0.03 -45.3 +/- 0.01 942.0 +/- 0.02
200 -51.4 +/- 0.03 5.0 +/- 0.03 -46.4 +/- 0.01 945.7 +/- 0.08
220 -50.7 +/- 0.03 5.5 +/- 0.03 -45.3 +/- 0.01 942.0 +/- 0.03
200 -53.9 +/- 0.05 5.0 +/- 0.03 -48.9 +/- 0.05 960.6 +/- 0.06
220 -53.4 +/- 0.06 5.5 +/- 0.03 -48.0 +/- 0.05 955.5 +/- 0.03
200 -54.0 +/- 0.05 5.0 +/- 0.03 -49.0 +/- 0.05 959.3 +/- 0.06
220 -53.5 +/- 0.06 5.5 +/- 0.03 -48.0 +/- 0.05 955.8 +/- 0.01
Pure Ice No E-Field
Pure Ice With E-Field
Ice + Ion with NO E-Field
Ice + Ion with E-Field

128

4.3.2 Structural Properties from the Radial Distribution Function
The RDF analysis of various ice lattices also validates excellent structural stability. The
structural comparison of RDF in lattices containing TIP4P, TIP4P/ice, TIP4P-Ew, and
TIP4P/2005 water models is shown on Figure 4.3.2.1. Only the TIP4P/ice model has the ideal
first peak separation distance at 0.275 nm (2.75 Å) while the first peak of other water models are
located at 0.270 nm. The location of the second peak of all four water models is the same at
about 0.450 nm. The first minima of all water models are at nearly 0 indicating excellent
separation between the first and second solvation shells as expected in solid phases. The
TIP4P/ice model lattice also has the highest first and second peak of 5.172 and 2.642,
respectively. Again, the ice lattice from the TIP4P/ice water model is the most tightly packed
structure further suggesting maximum stability relative to all TIP4P models.
The structural changes in g(O-O) of ice due to the addition of ion and electric field are
also analyzed. Due to slow ion diffusion, it is expected that the g(O-O) will be sensitive to initial
ion locations. This is substantiated by our simulations as the g(O-O) exhibits an appreciable
drop in packing efficiency with the additions of ion and electric field. As shown in Figure
4.3.2.2, the first peak of the ice structure with no ion and electric field is 5.14 as oppose to 4.72
for system with added ions and electric field. The first minima is also deeper for the pure ice
system indicating a more tighter bound first solvation shell with clearer separation from the
second solvation shell. The second peak of the pure ice system is also higher at 2.62 as opposed
to 2.14. Different ion configurations exhibited significantly less structural differences as shown
in Figure 4.3.2.3 after1 μ-sec simulations time. Ion configuration 1 exhibited the highest first
and second peaks of 5.23 and 2.24, respectively, but the other configurations are within 10% of
the peak values. The g(O-O) is also similar at temperatures 200K and 220K (see Figure 4.3.2.4).
Unexpectedly, the lower first minima indicated improved ice structure at 220K though the true
effect of temperature is most likely masked by the presence of ions. The radius of the first and
second solvation shell remain relatively constant at approximately 2.75 Å and 4.50 Å in all the
simulated systems, indicating the structure of ice is maintained throughout the simulation
regardless of the additions of ion and electric field.

129

Figure 4.3.2.1: The radial distribution functions of Oxygen-Oxygen comparing ice system
containing TIP4P, TIP4P/ice, TIP4P-Ew, and TIP4P/2005 water models.

Figure 4.3.2.2: Average radial distribution function of Oxygen-Oxygen comparing pure
ice system and ice + ion system. Temperature = 200K. The RDF of ice + ion with electric field
was calculated at the end of simulation (1 μ-sec).

130

Figure 4.3.2.3: Average radial distribution function of Oxygen-Oxygen comparing ice +
ion system with electric field using various initial ion configurations. Temperature = 200K. The
RDF was calculated after more than 1 μ-sec simulation duration.

Figure 4.3.2.4: Average radial distribution function of Oxygen-Oxygen comparing ice +
ion system with electric field at 200K and 220K simulation temperatures. The RDF was
calculated after more than 1 μ-sec simulation duration.

131

The RDF of O and H with the K
+
and Cl
-
ions revealed some stark differences in the local
structures near the ice oxygen and hydrogen atoms. After 1 μ-sec simulation duration as shown
in Figure 4.3.2.5a, the first peak g(O-K
+
) and g(O-Cl
-
) remain approximately 2.78 Å and 3.18 Å
apart, respectively. The difference is mostly attributable to the size differences between the K
+

and Cl- ions. Both first peaks are relatively sharp but the g(O-K
+
) is significantly higher than
that of g(O-Cl
-
) at 5.36 as opposed to 4.80 suggesting higher packing efficiency due to the
attraction between K
+
and partial negative charges on O. The second peak of both g(O-K
+
) and
g(O-Cl
-
) are located at about the same radius of approximately 4.70 Å. The g(H-K
+
) and g(H-Cl
-
) as shown in Figure 4.3.2.5b exhibited the first peak locations at 3.28 Å and 2.26 Å,
respectively. This is an indication that the attractive forces between the hydrogen and Cl
-
ions
are tighter than that between the oxygen and K
+
ions and is also attributable to a combination of
electrostatic forces and atomic size. The magnitude of the first peak of g(H-Cl
-
) is nearly twice
that of g(H-K
+
) at 4.5 versus 2.6. Combined with the observation that the first and second peaks
for g(H-K
+
) are also significantly broader than those of g(H-Cl
-
) suggests a significantly lower
packing efficiency of K
+
ions solvation shells around H as compared to the Cl
-
ions. The smaller
atomic size of H provides more available space for the preferential ion Cl- to approach closer,
and enhances the electrostatic attraction. Similar results were also obtained from the work of
Byrk and Haymet (2005). As a consequence of such observed lower packing order, the
likelihood of ionic motion within the ice lattice could increase for K
+
ions as compared to Cl
-

ions.

132

Figure 4.3.2.5a: Average radial distribution function of Oxygen-K
+
ion and Oxygen-Cl
-
ion.
The RDF was calculated after more than 1 μ-sec simulation duration.

Figure 4.3.2.5b: Average radial distribution function of Hydrogen-K
+
ion and Hydrogen-Cl
-
ion.
The RDF was calculated after more than 1 μ-sec simulation duration.

133

From the g(K-Cl), g(K-K), g(Cl-Cl), we also attempted to analyze the impact of the
electric field in the ion structure throughout the simulation. As shown in Figure 4.3.2.6, the
initial g(K-Cl) exhibited the first peak of 70 indicating high packing efficiency from strong
electrostatic forces. The addition of electric field tends to weaken the electrostatic forces as
observed by the decrease the first peak by more than 50% down to 30.0 after 1 μ-sec. This is
physically accurate as ion tends to move apart due to the added force from the electric field. A
relatively constant 3.0 Å location of the first peak throughout the entire simulation time suggest
no ion-paring in the ice structure. The g(K-K) and g(Cl-Cl) indicate significant structure
changes with the addition of the electric field. For K
+
ion shown in Figure 4.3.2.7, the electric
field decrease the location of the first peak by approximately 0.3 Å and the second peak by
approximately 0.5 Å. The magnitude of the second peak nearly double under the influence of the
electric field indicating like ions are moving along the same direction. Figure 4.3.2.8 for g(Cl-
Cl) also exhibited small peak formations of less than 1.0 occurring at the separation distance of
5.5 Å which suggests the onset structure changes as well. The magnitude of the second peak at
6.7 Å dramatically increased with simulation time under the influence of the electric field and
suggests secondary structural formation as well. Even after 1 μ-sec simulation time, such a
dynamic process does not appear to have reached equilibrium under the influence of the electric
field.

Figure 4.3.2.6: Average radial distribution function of K-Cl calculated at various simulation
times. The radial distribution functions with the E-field are also time-averaged over 50 ns
simulation time.

134

Figure 4.3.2.7: Average radial distribution function of K-K calculated at various
simulation times. The radial distribution functions with the E-field are also time-
averaged over 50 ns simulation time.

Figure 4.3.2.8: Average radial distribution function of Cl-Cl calculated at various
simulation times. The radial distribution functions with the E-field are also time-
averaged over 50 ns simulation time.

135

4.3.3 Ionic Transport Mechanism and Jump Events
Details of the ion transport mechanism are initially observed through individual ion
displacement profiles. The displacement profile of each individual ion throughout the simulation
was calculated in reference to the initial position at t = 0. Figure 4.3.3.1A,B summarized the
dynamic displacement profiles of a representative set of K
+
and Cl
-
ions; it is observed that
numerous ions experience sudden large discreet changes in displacement which can be described
as “jump” events. With the understanding that the ice molecules are mainly stationary, such
ionic jump event must occur across the ice cavities. The ions also exhibit period of little or no
displacement as described by the long plateau; this describes “caged” motion. At a glance, the
K
+
ions appear to exhibits higher overall displacement throughout the simulation time as well as
higher frequency of jump events than that of the Cl
-
ions.
Meaningful statistics and dynamics of the jump events across the ice lattice cavities are
characterized using a similar technique developed for other crystalline materials.
83
This
methodology entails the calculation of the running average and its corresponding standard
deviation of the displacement over small simulation time window (<10000 ps). As illustrated by
the example comparison of the running average versus actual displacement profile shown in
Figure 4.3.3.2a, the moving average provides smoothing of the raw data thereby enabling
detection and counting of jump events. Jump events are defined by the criteria of displacement
within a small time-step (<15000 ps) exceeding twice the standard deviation over the entire
simulation (2*σ), which is a self-consistent value determined from the mean-square displacement
in continuous random-walk processes. An example is illustrated in Figure 4.3.3.2b. This
technique not only provides critical information on individual ion jump frequencies over long
simulation times, but also information on the jump dynamics such as jump timescale (defined as
the time between successive jumps).

136

Figure 4.3.3.1: Representative displacement profiles for K+ ions throughout the
simulations.

Figure 4.3.3.1: Displacement profiles for throughout the simulations: (A) representative set
of K+ ions, (B) representative set of Cl- ions.

(B)
(A)

137

Figure 4.3.3.2: Methodology of ionic jump characterization: (A) comparison of actual vs
running average of displacement profiles, (B) running average of displacement with
corresponding standard deviation profiles. The “jump” threshold = 2*σ defines a jump event
(red circles).

(A)
(B)

138

To acquire meaningful jump statistics, the accumulative jump count for individual K
+
and
Cl
-
ions are determined over production simulation runs in the presence of the electric field. The
total frequency of jump events from all three ion configurations are compared in Figure
4.3.3.3A,B. At 200K, the jump count of both K
+
and Cl
-
ions exhibited excellent repeatability
from the 3 initial ion configurations, suggesting that jump frequencies are relatively independent
initial ion locations attributable to the slow diffusion processes in solid-state materials. Table
4.3.3.1 summarizes the jump rates K
+
and Cl
-
ions over the simulation time range of 200000 ps
to 1 μ-sec. The average jump rate for K
+
is found to be approximately 70 jumps per μ-sec, a
factor of 5 times increase over the jump rate of Cl- ions (15 jumps per μ-sec). The jump rates
were also calculated at 220K for only the initial ion configuration 2 system. Relative to 200K,
the jump rates for both K
+
and Cl
-
ions increase similarly at 220K by a factor of 3.2 and 4.8,
respectively. Therefore, the presence of ions, most notably K
+
ions, and added kinetic energy
appear to contribute substantially to the transport properties of ice as observed by jump events.
The time-scale of the jump events was also evaluated in terms of ionic type and simulations
temperatures. At 200K, the average jump time-scale over the entire simulation for K
+
and Cl
-

ions are 163.8 ns and 235.3 ns, respectively. At 220K, the average jump time-scale over the
entire simulation for K
+
and Cl
-
ions decreased to 81.3 ns and 115.1 ns, respectively. Once
again, the increases in temperature plays a critical role in jump frequency but the distribution of
time-scales vary over a wide range suggesting random jump events occurrences. Figure
4.3.3.4A,B summarizes the probability distribution functions of the jump time-scale for both K
+

and Cl
-
ions at both 200K and 220K from a 16 ns time average window. The highest jump
probability resides in the jump time-scale of 100000 ps or less, and significantly decreases with
increasing jump time-scales. Again, the highest probability of jump time-scale also increases
with increasing temperature by 0.20 or higher. Similar results were obtained from the
distribution jump time-scales of 100 ns or less using a substantially smaller time average window
of 0.5 ns in order to capture jump events occurring within smaller time steps. As shown in
Figure 4.3.3.5, the highest probability of jump events also occurs at the lowest time-scale. This
explicably suggests that the local ice structure surrounding ions plays an important role, and is
most likely attributable to the presence of water-like structures facilitating free ions motion.

139

Figure 4.3.3.3: Total jump counts from simulations at T = 200K of 3 initial ion
configuration. (A) K
+
ion jump count, (B) Cl
-
ion jump count.

(A)
)

(B)
(B)

140

Table 4.3.3.1: Summary of Jump Rates of all 3 ion configurations at T=200K and 220K.

Jump Rate (Jumps/ μ-sec)
& Jump Rate Ratios
Ion
Configuration 1
Ion
Configuration 2
Ion
Configuration 3
Average
K
+
(200K)
61.9 67.4 81.1 70.1
Cl
-
(200K)
14.0 18.7 12.1 14.9
K
+
(220K)
-- 214.5 -- --
Cl
-
(220K)
-- 90.5 -- --
Ratio Jump Rate
K+/Cl-
4.4 3.6 6.7 4.9
Ratio Jump Rate K
+
220K/200K
-- 3.2 -- --
Ratio Jump Rate Cl
-
220K/200K
-- 4.8 -- --
(A)
)

(B)

141

Figure 4.3.3.4: Jump time-scale distribution function at different temperatures over 1 μ-sec
simulation time: (A) K
+
ion jump time-scale distribution function, (B) Cl
-
ion jump time-scale
distribution function. Running average time step of 16 ns.

Figure 4.3.3.5: Jump time-scale distribution function K
+
and Cl
-
ions over 100000 ps
simulation time at T=200K. Running average time step of 0.5 ns.

(B)
)

(B)

142

In order to probe the ion interaction energy between the water molecules and the ions, we
performed a local energy analysis within the ice lattice. Such energy profiles provide
information on the energy barriers experienced by the ions across the ice lattice in an effort to
understand its transport mechanisms and potential migration pathways. The total pair potential
equation is described by
   
   











 
N
i
N
i j ij r
j i
N
i
N
i j ij
ij
ij
ij
r
q q
r
B
r
A
E
1 1
6 12
) ( 332

4.3.3.1
where N is the total number of charged particles containing center, r
ij
is the separation radius
between particle i and j. The parameters A and B in equation 4.3.3.1 are defined as



6
6
12
2
2



o
o ij
o ij
r
r B
r A

where ε
r
is the dielectric constant of 1, the value of 332 is a unit conversion factor, and the ε and
σ are standard Lennard-Jones parameters describing the depth of potential well and finite
distance of zero inter-particle potential, respectively. For these calculations, we employed the
same force-field parameters corresponding the TIP4P water model as AMBER.
84,85,86
Both local
potential energy and energy contour calculations employ the same parameters a summarized in
Table 4.3.3.1.
Table 4.3.3.1: Summary of parameters for pairwise energy calculation.

M = Charge point of TIP4P model, K = Potassium ion, Cl = Chloride ion, H = Hydrogen
Pair Interaction or
Charge
A
parameter
B
parameter
q
M-M 5.999E5 6.099E2 --
K-M 4.050E5 1.075E2 --
Cl-M 4.202E6 1.447E3 --
q(M) -- -- +1.04
q(H) -- -- -0.52
q(K) -- -- +1.0
q(Cl) -- -- -1.0

143

Using the trajectory coordinates throughout the simulation duration, the local potential
energy with the water molecules was determined within a 1.6 nm cut-off radius from each ion.
A representative local energy profiles of K
+
and Cl
-
ions containing jump events were plotted
against their corresponding displacement profiles as shown in Figures 4.3.3.6A,B and
4.3.3.7A,B, respectively. The reported local potential energies are normalized with the total
number of water within the cut-off radius. The range of local potential energy is approximately
+/-130-134 KJ/mole which is consistent with other studies (Kadota et al, 2007). As illustrated in
these plots, the ion jump events have a direct correlation to its local potential energy profile. It is
observed that the decrease in local energy of both K
+
and Cl
-
ions directly corresponds to the
occurrences of the jump events. However, no conclusions can be made on whether the local
energy barriers initiate such jump events. It is also observed that the magnitude of the local
energy changes before and after the jump events are significantly higher at 220K as compared to
200K, thus, kinetic energy does play a significant role on jump events. Between energy profiles
of K
+
and Cl
-
ions, there exist no appreciate differences in the magnitude of the local energy
changes before and after jump events. Therefore, it is the number of local water molecules
surrounding the ions that play a critical role in ion migration pathways.
Using Equation 4.3.3.1, a grid-based interaction energy calculation technique was also
performed to generate energy contour profiles of K
+
and Cl
-
ions. The energy contour algorithm
computes the total energy at each of the 10
6
grid positions (0.1 Å apart) within the cubic box
length of 10 Å from arbitrary initial coordinates of an ion within the ice lattice. The LJ energy
summation was computed between O atom positions within a cut-off distance of 40 Å. To
account for long range effect, the electrostatic interactions employ a cut-off distance of 15 Å
using standard Ewald summation method for rigid molecules. An idealized ice lattice with no
energy minimization containing 1600 water molecules was utilized to generate these energy
contour profiles. To ensure the energy calculations are consistent over a global scale in the
lattice, numerous energy maps over a 10 Å cubes of ice are performed at different locations of
the lattice.

144

Similar conclusions could be drawn in the contour maps of interaction energies between
the ion and the water molecules. A representative cross-sectional energy contour plots are
shown in Figure 4.3.3.8A,B. The energy calculations were performed by moving the K
+
and Cl
-

ions incrementally by 0.2 Å across same 10 Å x 10 Å x 10 Å ice lattice space. Though energy
minima for both ions are located in the ice cavities, the K
+
ion exhibits a lower intermolecular
energy gap between adjacent ice molecules than that of the Cl
-
ion. This observation is also
attributable to ionic size as the Cl
-
is significantly larger than the K
+
ions. In addition, both
contour slices exhibited the same behavior suggesting a global characteristic over entire ice
lattice. This lower interaction energy of the K
+
ions supports the notion of wider migration
pathways which could also facilitate jump events across the adjacent ice cavities. Therefore, this
offers an additional explanation for the observed higher frequency of jump events as well as
longer displacement distances over the same simulation duration by the K
+
ions as compared to
the Cl
-
ions.

145

Figure 4.3.3.6: Comparison of displacement running average and calculated local
interaction energies for K+ ions: (A) T = 200K, (B) T = 220K. Local energy calculations are
based on 6000 ps running average. Negative energy values are due to positive charge K
+
atoms.

(B)
)

(B)
(A)
)

(B)

146

Figure 4.3.3.7: Comparison of displacement running average and calculated local
interaction energies for Cl
-
ions: (A) T = 200K, (B) T = 220K. Local energy calculations are
based on 6000 ps running average. Positive energy values are due to negative charge Cl
-
atoms.

(A)
)

(B)
(B)
)

(B)

147

Figure 4.3.3.8: Cross-section of interaction energy contour maps: (A) Z = 17.757 Å, (B) Z
= 25.957 Å. The computation is performed over 10 Å box with a 0.2 Å step. The cut-off
distances for the Lennard-Jones and electrostatic forces are 40 Å and 15 Å, respectively. Large
white spots surrounded by red circle are locations of infinite energy and should be colored red.
Both energy contour maps suggest the K
+
ions have a wider migration channels than Cl
-
ions
over the same ice lattice space.
Lx (Angstroms)
Ly (Angstroms)

16 18 20 22 24
17
18
19
20
21
22
23
24
25
0
1
2
3
4
5
6
7
8
9
x 10
5
Z = 17.757 A
Lx (Angstroms)
Ly (Angstroms)

16 18 20 22 24
17
18
19
20
21
22
23
24
25
0
2
4
6
8
10
12
14
16
18
x 10
5
Z = 17.757 A
Lx (Angstroms)
Ly (Angstroms)

20 22 24 26 28
22
23
24
25
26
27
28
29
30
0
1
2
3
4
5
6
7
8
9
x 10
5
Z = 25.957 A
Lx (Angstroms)
Ly (Angstroms)

20 22 24 26 28
22
23
24
25
26
27
28
29
30
0
2
4
6
8
10
12
14
16
18
x 10
5
Z = 25.957 A
(B) Energy contour profiles for K
+
(left) and Cl
-
(right) over
same the ice lattice space at Z = 25.957Å

(A) Energy contour profiles for K
+
(left) and Cl
-
(right) over same
the ice lattice space at Z = 17.757Å

148

4.3.4 Ionic Transport Properties and the Mean-Square Displacement
The contribution of ion transport in ice is quantified from MSD analysis. The equation for
MSD is defined by
2
) 0 ( ) ( r t r MSD  
4.3.4.1
where r(t) is the coordinates of a particle at time t, r(0) is the coordinates of the particles at initial
time 0, and <..> denotes the averaging over all the particles of interest. We calculated the
average MSD of water molecules and ions in order to examine their rates of diffusion throughout
the simulation. The MSD profiles highlighting the effect of electric field on pure ice simulations
are shown in Figures 4.3.4.1 and 4.3.4.2. In pure ice, the MSD is significantly higher despite
shorter simulation duration than the simulation without the electric field. This suggests ice
mobility noticeably decreases under the influence of the electric field which is similar to the
behavior observed in liquid water (Wei et al., 2005). The same is also observed with simulations
comparing the ice doped with ions where the MSD also decreased under the influence of the
electric field. It is also observed that ice doped with ions exhibits higher rate of diffusion than
pure ice with or without the applied electric field. When comparing the average MSD of K
+
and
Cl
-
ions versus only the ice molecules in the same simulation system, K
+
ions exhibited a highest
MSD as shown in Figure 4.3.4.3. The average MSD of pure ice and Cl
-
ions are nearly identical
with not electric field present. Under the influence of the electric field, the average MSD of ice
is lower than both the MSD of K
+
and Cl
-
ions. The influence of the electric field decreased the
average MSD of all ice and K
+
and Cl
-
ions as compared to the MSD with no electric field
present. Therefore, MSD results substantiates the effect of the electric field on ion diffusion and
the dominate role of Cl- and especially K
+
ion motion on the transport properties of ice.
The average MSD of only the K
+
and Cl
-
were calculated along the individual X,Y,Z-axis to
provide further understanding on the translational motion caused by the applied electric field.
These calculations were performed by taking only the X,Y,Z-components of the particle
coordinates of Equation 4.3.4.1. At 200K, the electric field applied along the Z-axis significantly
increased the average K
+
and Cl
-
MSD along the same direction as shown in Figure 4.3.4.4 as
compared to the X and Y-axis. This illustrated by the slope of the MSD along the Z-axis which
exhibits an higher rate of increase throughout the simulation in comparison with the MSD slope

149

along the X,Y-axis. At 220K, however, not appreciable impact was observed in the average
MSD along the X,Y,Z-axis at the end of the simulation as shown in Figure 4.3.4.5. Though the
MSD along the X-direction is higher, the slope of the MSD along the direction of the electric
field appears higher after 400000 ps simulation time. This suggests that the added kinetic energy
at 220K is an important factor in masking the influence of the relatively weak applied electric (1
1 V/cm). Residence time maybe also a consideration on the influences of an electric field as it
may take more time for the ions to polarize in the direction of the applied field to overcome the
added kinetic energy from higher temperatures. Overall, our results indicate that temperature
and kinetic energy remains the dominate factor in ion transport mechanisms in the ice lattice. It
is presumed that ion jump events are facilitated with increases in temperature with the migration
path dictated by local ice structural changes as previously discussed. Therefore, the additional
driving force provided by the applied electric field play appears to play lessor role in ion
transport at higher temperatures as compared to lower temperatures.

150

Figure 4.3.4.1: Average MSD of pure ice with and without electric field. Simulation
temperature = 200K. Simulation time is normalized to 1 for comparison purposes. Simulation
durations with electric field are typically more than 10 times longer than the simulations without
electric field.

Figure 4.3.4.2: Average MSD of ice and ion with and without electric field. Simulation
temperature = 200K. Simulation time is normalized to 1 for comparison purposes. Simulation
durations with electric field are typically more than 10 times longer than the simulations without
electric field.

151

Figure 4.3.4.3: Average MSD comparing K
+
ions, Cl
-
ions, and water in the same
simulation lattice: (Left) without the applied electric field, (Right) with the applied electric field.
Ion configuration lattice 2. The simulation with electric field is Simulation temperature = 200K.

Figure 4.3.4.4: Average MSD of K
+
(Left) and Cl
-
(Right) ions with the applied electric
field along the X,Y,Z-axis. Ion configuration lattice 2. Simulation temperature = 200K.

152

Figure 4.3.4.5: Average MSD of K
+
(Left) and Cl
-
(Right) ions with the applied electric
field along the X,Y,Z-axis. Ion configuration lattice 2. Simulation temperature = 220K.

153

The translational diffusion coefficient, D, can be obtained from the integral of the velocity
autocorrelation function or from the MSD using the expression

2
) 0 ( ) ( lim
6
1
i i
t
r t r
dt
d
D  
 
4.3.4.1
where r
i
is the positional vector. A comparison of the self-diffusion coefficient (D
s
) of pure ice
from our simulation and experimental studies is plotted in Figure 4.3.4.6. The experimental
results at 160K and 263.15K come from the work of Brown and George (1996) and Ramseier
(1967), respectively, with values of D
s
directly proportional to temperature. From MD
simulations with the addition of the electric field, the value of D
s
at T = 200K is approximately
6.30 x 10
-14
cm
2
/sec and increases with temperature to approximately 8.1 x 10
-13
cm
2
/sec at
220K. Overall, our D
s
values from simulation compares well with that of experimental findings.
Though it is not strictly defined as diffusion due to the polarizing influence of the electric
field, we used the diffusion coefficient to characterize the effect of ion motion in simulations
containing the electric field. The diffusion coefficients from simulations of K
+
and Cl
-
ions in an
electric field are summarized in Table 4.3.4.1. The slope of the diffusion coefficients are taken
at the simulation time between 400000 ps to 800000 ps from the average MSD. At 200K, the
MD simulation from all three initial ion configuration exhibited similar ice diffusion coefficients
with an average of 1.91 x 10
-10
cm
2
/sec. Despite having a substantially higher initial ion
separation radius of greater 10 Å, initial ion configuration 2 also has the highest diffusion
coefficient of ice at 2.980 x 10
-10
cm
2
/sec. This suggests that no deliberate migration pathway
was created during initial ion substitution insertion process. In addition, the average diffusion
coefficient of all three initial configurations with the applied electric is about 1 order of
magnitude lower than those without the influence of the electric field at 5.442 x 10
-10
cm
2
/sec.
From the estimated diffusion coefficient of the individual ions, it is determined that the average
K
+
diffusion coefficient is nearly one order of magnitude higher as compared to that of the ice or
the Cl
-
ions. The average diffusion coefficient of K
+
ions is 1.450 x 10
-9
cm
2
/sec as oppose to
4.66 x 10
-10
cm
2
/sec for Cl- ions under the influence of the electric field. Similar results were
observed from the T = 220K simulation. Overall, the diffusion coefficient suggests that ion
migration is the main contributor in ice lattice. While pure ice does exhibit some transport most

154

likely due to defect motion, it is more likely the ions (especially K
+
ions) that dominate such
transport properties even under the influence of the electric field as support by the observed
cavity jump events and their corresponding diffusion coefficients.
From the estimated diffusion coefficients along the X,Y,Z-axis, we also determined the
effects of the electric field and temperature on transport mechanisms in ice. These calculations
were also taken at the simulation time between 400000 ps to 800000 ps from the average MSD.
At the lower temperature of 200K, the applied electric exhibits a stronger influence on the
diffusion coefficient. The diffusion coefficient for both K
+
and Cl
-
ions is about 1 order of
magnitude higher along the direction of the applied field (Z-axis). As displayed in Table 4.3.4.2,
the Dz (diffusion coefficient along the Z-axis) is 1.236 x 10
-9
cm
2
/sec for the K
+
ions and 9.487 x
10
-10
cm
2
/sec for Cl
-
ions. At 220K, the diffusion coefficients along the X,Y,Z-axis is nearly one
order of magnitude (factor of 10) higher than those at 200K. However, the Dz value along the
direction of the applied electric field is only higher by a factor of 2 at 4.025 x 10
-9
cm
2
/sec the K
+

ions. For the Cl
-
ions, no appreciable polarization due to the electric field is observe as both Dx
and Dz are nearly identical at 1.699 x 10
-9
cm
2
/sec and 1.965 x 10
-9
cm
2
/sec, respectively. As
indicated by these results, higher temperature appears to mask the influence of the electric field
causing ions to migrate more randomly as oppose to the direction of the applied electric field
observed at the lower temperature of 200K.

155

Figure 4.3.4.6: Comparison of pure ice self-diffusion coefficients from simulations and
experiment. The Ds values at 160
o
K (Brown and George, 1996)
87
and 263.15
o
K (Ramseier,
1967)
88
provide comparison basis for our simulation values.

Table 4.3.4.1: Summary of diffusion coefficients from simulation with electric field on all
three ion configurations at T = 200K.

*Diffusion coefficient of ice is based on avg MSD of ice, K
+
, and C
+
ions.

Table 4.3.4.2: Summary of K
+
and Cl
-
ion diffusion coefficients along the X,Y,Z-axis from
simulations with electric field in Z-direction.

Type
Ion Config 1
(with E-Field)
Ion Config 2
(No E-Field)
Ion Config 2
(with E-Field)
Ion Config 3
(with E-Field)
Average
(with E-Field)
Stan Dev
(with E-Field)
Ice 9.940E-11 5.442E-09 2.980E-10 1.750E-10 1.908E-10 1.002E-10
K
+
ions only
1.738E-09 1.060E-08 1.745E-09 8.605E-10
1.448E-09 5.087E-10
Cl
-
ions only 1.689E-10 4.604E-09 9.468E-10 2.831E-10
4.662E-10 4.200E-10
Diffusion Coefficients (cm
2
/sec)
Temperature (K) D
X
Dy Dz D
X
Dy Dz
200 3.346E-10 2.289E-10 1.236E-09 1.534E-10 8.284E-11 9.487E-10
220 2.142E-09 1.172E-09 4.025E-09 1.699E-09 5.587E-10 1.965E-09
K
+
Ion Diffusion Coefficients (cm
2
/sec) Cl
-
Ion Diffusion Coefficients (cm
2
/sec)

156

4.3.5 Electrical Properties
The electrical property of relative permittivity (or dielectric constant) was determined
from simulations to offer direct comparison with experimental infinite dielectric constant, ε
∞
.
The dielectric constant, ε
r
, from simulation of ice is determined by the magnitude and direction
of the molecular dipole moments. From the Ewald approach, ε
r
is defined by






  
2
2
3
4
1 M M
T k V
o b
r
 



4.3.5.1
 
2
2
2
2 2 2
2
z y x z y x
M M M M M M
M M M M M
    
    

where ε
o
is infinite dielectric constant of 1, k
b
is the Boltzmann constant, V is the average
volume, T is the average system temperature, and
2
M
-
2
M
is the fluctuation of the total
system dipole moment as calculated using Equation 4.3.5.2. The profiles of the dielectric
constant of pure ice at 200K and 220K from simulation are shown in Figure 4.3.5.1. After only
approximately 50000 ps simulation time, the dielectric constant has equilibrated even under the
influence of the applied electric field. The equilibrium value of the dielectric constant at 200K
(1.36) is very similar to that at 220K (1.39) which indicates relatively temperature independent.
The overall observed low values of dielectric constant appears to be influenced by the presence
of electric field in a form of higher restrain rotational motion of individual water molecules along
the direction of the applied field. The effect of the electric field on systems of ice with ions is
shown in Figure 4.3.5.2. As expected, the dielectric constant is higher without the influence of
the electric field from more random motion of water molecules. With the addition of ions, the
dielectric constant is significantly higher than of pure ice. As shown in Figure 4.3.5.3, all three
ion configurations at 200K exhibited similar dielectric constant of nearly 2.0 or higher. At 220K,
the dielectric constant has increased to 4.0 and beyond after 400000 ps simulation time. This is
in contrary to what is expected but is consistent with the observed increases in jump frequency
and diffusion coefficients at higher temperatures which tend to facilitate molecular motion due to
local melting near the ions.

157

Figure 4.3.5.1: Comparison of dielectric constant profiles from simulations of pure ice systems
with applied electric field at 200K and 220K. The equilibrium value of the dielectric constant is
1.36 and 1.39 for 200K and 220K, respectively.

Figure 4.3.5.2: Comparison of dielectric constant profiles from simulations of ice and ion
systems with and without the influence of the electric field at 200K.

158

Figure 4.3.5.3: Comparison of dielectric constant profiles from simulations of ice and ion
systems with applied electric field at 200K and 220K over 1 μ-sec simulation time.

The direct comparison of experimental and simulation dielectric constant is plotted in
Figure 4.3.5.4. As shown here, the simulated dielectric constant is relatively close to that of the
experimental values. In pure ice, the dielectric constant from simulation is relatively
independent of temperature at 1.35 compared to the experimental dielectric constant of 3.0. In
ice simulations with ions, however, the dielectric constant is shown to be temperature dependent.
The values of 2.39 and 4.38 were determined after 1 μ-sec simulation time at 200K and 220K,
respectively. Large disparity of values can be attributable to the differences in the experimental
conditions and ice nucleation processes. The ion distribution from simulation is substantially
more controlled as compared to experiments where ions tend to move away from the sense
probes during the freezing process (see Section 3.3). This also substantiates the similarity of
experimental dielectric constants at different ionic concentrations. Extrapolating out to a lower
temperature of 200K, the dielectric constant is substantially lower than that from experiment.
Though experimental values of dielectric constant show conflicting results in differences with
thermal energy, the dielectric constant from simulation represents a more ideal behavior of
dielectric constant of ice doped with ions due to minimal experimental factor.

159

Figure 4.3.5.4: Comparison of experimental versus simulation dielectric constant. The
experimental dielectric constants are extracted at 100 kHz. The simulated dielectric
constant is extracted at the end of the simulation. The value at T = 200K is the average of
three simulations from the 3 ion configuration.

160

As a final direct comparison, the conductivity of ice doped with ions is calculated using a
generalized form of the Nernst-Einstein equation

i i
i
i
b
D q n
T k
e
2
2

 
4.3.5.2
where n
i
is the number of i ions in the simulation system, q
i
is the charge of the ion, D
i
is the
diffusion coefficient, and e is the elementary charge. For purposes of simplifying the
calculations, the diffusion coefficient was determined by the averaged MSD over all the same
ionic species. Figure 4.3.5.5 compares the experimental and simulated conductivity of the ice
doped with ion. Similar to experimental conditions, the simulation systems are also under the
influence of the electric field in conductivity measurement. Here, the simulated conductivity
includes only the contributions from the ions since the contribution from ice is negligible. At
220K, the simulated conductivity is 3.4 x 10
-7
S/cm as compared to the experimental value of
approximately 5.0 x 10
-8
S/cm. Extrapolating out to a lower temperature of 200K, the simulated
conductivity is 4.4 x 10
-8
S/cm averaged over all three ion configurations. Therefore, both
simulated and experimental conductivities exhibit a decrease with decreasing temperature which
is consistent with physical expectations.

161

Figure 4.3.5.5: Comparison of experimental versus simulation conductivity. The simulated
conductivity was calculated with the ions. The value at T = 200K is the average of three
simulations from the 3 ion configuration.

162

4.3.6 Visualization and Trajectory Plots

Finally, we employed numerous variations of trajectory plots to provide detailed
examination and visualization to support the mechanisms of ion transport in the ice lattice. First,
trajectory snapshots at the start and at the end of simulation time are provided to assess the
integrity of ice lattice as well as to provide some observations of the local environment around
individual ions. The electric field was added after t = 0 ns. The large circular dots represent the
K
+
and Cl
-
ions within the lattice. At T = 200K and 220K, it is shown that the hexagonal ice
lattice structure is retained after 1000 ns or more simulation time under the influence of the
electric field (Figure 4.3.6.1 and Figure 4.3.6.2). The ions at 200K do not appear to migrate
appreciably after more than 1000 ns simulation time as shown by the 3D plot. At 220K,
however, a significant number of ions do appeared to have migrated substantially from their
initial positions after 1000 ns simulation time. A close-up examination of local environment of
surrounding water molecules around the ions shown in Figure 4.3.6.3 also exhibits a more
disordered structure compared to lattice locations with ions. This higher quantity of dislocated
water molecules from the ice lattice positions are defects and can ultimately be interpreted as
protonic-types as proposed by previous investigators (Bjerrum, 1952). The distribution of
displaced water molecules extend across many ice cavities which suggests a probable cause is
due to the migration pathway of the individual ions.

163

Figure 4.3.6.1: Trajectory 3D snapshot at t = 0 simulation time and t = 1040 ns simulation
time at T=200K. The electric field was applied after t = 0 simulation time. The green and
red lines are the water molecules. The large white and blue dots are the centers of K
+
and
Cl
-
ions.

Figure 4.3.6.2: Trajectory 3D snapshot at t = 0 simulation time and t = 1000 ns simulation
time at T=220K. The electric field was applied after t = 0 simulation time. The green and
red lines are the water molecules. The large white and blue dots are the centers of K
+
and
Cl
-
ions.

164

Figure 4.3.6.3: Trajectory snapshot after t = 300 simulation time at T=220K. (Right) 2D
ice lattice with close-up portion in the black square shown on the left picture. (Left) The
blue circle highlights highly disordered water structures surrounding the ions.

The total trajectory pathway of all individual ions at T = 200K and 220K are shown in
Figures 4.3.6.4 and 4.3.6.5. By projecting the trajectories over both XY-plane and YZ-plane, we
were able to obtain a more accurate depiction of ion transport mechanism within the ice lattice
throughout the simulation time. Due the large lattice dimension of 53.9 Å x 62.2 Å x 58.7 Å, it
difficult to illustrate significant translational motion at T = 200K. The majority of ions (both K
+

and Cl
-
) appears stationary or trapped at the relatively same vicinity within its initial insertion
positions of the ice cavities. Ion migration is more clearly illustrated at the 220K temperatures.
Over the same XY-plane and YZ-plane projections, the ions exhibit significantly higher
migration distances as indicated by a wider footprint of the trajectories throughout the simulation
as compared to simulation at T = 200K. Similarly, many ions during the 220K simulation
appear to remain within the same ice cavity as well. We also noticed that the migration path
even within a single cavity of all the ions appears to traverse along same direction which can be
attributable to the direction and the influence of the electric field.
Individual K
+
and Cl
-
ion total trajectories were plotted to further substantiate the presence
of ion migration pathways across multiple ice cavities. As shown in Figure 4.3.6.6 plot of the
total trajectories of a K
+
ion, jump events across multiple ice cavities at 220K temperatures is
detected over such a long simulation time. Similarly, the Cl
-
ions also jump across multiple ice
cavities as shown in Figure 4.3.6.7. From these trajectory plots, both K
+
and Cl
-
ions exhibited

165

multiple jump events across ice cavities. These jumps appear to take a substantial amount of
time to initiate as a majority of trajectory positions are lumped near the same positions as
illustrated by the large clump of overlaying trajectories. These bottleneck positions prior to such
jump events appear to reside either in ice cavities or along the ice lattice position. From these
observations, it is more apparent that ions incorporated in the lattice are the main facilitators of
defect formation as observed by the local disordered ice structures near the ions. More
importantly, the ions themselves also exhibit translational motion across multiple ice cavities
which can significantly contribute to the transport properties of ice as previously described.

166

Figure 4.3.6.4: Trajectory plot of all 60 ions over entire simulation time t = 0 to 1000 ns at
T=200K. (a) Projected on XY-plane. (b) Projected on YZ-plane. The Green spots are the
total trajectories of the K
+
ions. The Grey spots are the total trajectories of the Cl
-
ions.

(a)
(b)

167

Figure 4.3.6.5: Trajectory plot of all 60 ions over entire simulation time t = 0 to 1000 ns at
T=220K. (a) Projected on XY-plane. (b) Projected on YZ-plane. The Green lines are the
total trajectories of the K
+
ions. The Grey lines are the total trajectories of the Cl
-
ions.

(a)
(b)

168

Figure 4.3.6.6: Trajectory plot projected of XY (Left) and YZ (Right) planes of single K
+

ion over entire simulation time t = 0 to 1000 ns at T=220K. The blue lines are the ice
lattice. The green line is the center of trajectories on the K
+
ion.

Figure 4.3.6.7: Trajectory plot projected of XY (Left) and YZ (Right) planes of single Cl
-

ion over entire simulation time t = 0 to 1040 ns at T=220K. The purple are makeup the ice
lattice. The red line is the center of trajectories on the Cl
-
ion.

169

4.4 Conclusions
In an effort to further understand the physical and chemical processes emanating from the
molecular level, molecular dynamics simulations were performed to compliment the
experimental measurements on the electrical properties of ice. Prior to large-scale lattice
development of ice, we initially performed ice simulations comprising of only 1600 water
molecules of various types of potential water force-fields (TIP4P, TIP4P-Ew, TIP4P/Ice, and
TIP4P/2005) in order to determine the most suitable water model for ice simulations. We
determined that all water models are suitable for ice simulations. Though TIP4P/Ice lattice
exhibited the lowest total energy, the closest total energy to that reported in literature is TIP4P at
-46.4 kJ/mole. In terms of other physical properties, the heat capacity and density of the TIP4P
lattice are determined to be about 50.0 J/K-mole and 946 kg/m
3
, respectively, and is in good
agreement with experimental studies as well. Consequently, our larger simulation ice lattices
incorporate the TIP4P water model comprising of 6144 water molecules. Numerous types of
atomistic simulation were also employed in order to mimic the actual experimental conditions,
including the baseline pure ice system with and without an electric field, and KCl-doped ice with
and without an electric field.
Significant molecular details on the mechanisms of ice transport were observed in our
MD simulations. Ions play a substantial role in ice physical properties as potential energy
fluctuations increased by a factor of 2 in simulations of ice and ions over pure ice systems. The
density of ice increased with the addition of ion which also suggests significant structural
dynamics. From RDF analysis, we found significant differences between the local structures of
K
+
and Cl
-
ions with the ice molecules. The Cl
-
ions are packed closer to the water molecules
than the K
+
ions which is attributable to smaller atomic size and tighter electrostatic forces
between opposite charges. In addition, the RDF of K-Cl under the influence of the electric field
also initially increases the packing efficiency as observed by the first peak but subsequently
decreases with increasing simulation time which suggests a high degree of structural dynamics
surround the ions throughout the simulation. The RDF of K-K and Cl-Cl also indicated
significant structural changes throughout the simulation under the influence of the electric field.
From ion displacement profiles, we were able to describe the mechanisms of ion transport
in ice through jump events across adjacent ice cavities. The quantification of ionic jump events

170

and jump frequency offered great details on the transport mechanism especially the dominance
of K
+
ions over Cl
-
ions. The jump rate of K+ ions is a factor of 4 or higher than the Cl- ions and
is consistent in 3 separate simulations containing a unique initial ion configuration. More
interestingly, the distribution function of jump timescale of both ion types is highest at a short
timescale of 10000 ps or less. This strongly suggests that the jump events occur over short
timescales once it initiates, and hence it is strongly attributable to the local structure surrounding
the ions. From the trajectory plots, we were able to confirm ions “hopping” across ice cages
albeit on only selective individual ions (not on all ions).
The analysis of local interactions energies also supports such a mechanism for ion
transport. Though it is inconclusive on whether local potential energy gradient initiates the
observed jump event, we found such jump events directly correspond to the decrease in local
potential energy. The difference in the local potential energy between K
+
and Cl
-
ions is
negligible but the magnitude of the energy gradient during the jump events increases with
temperature due to the added kinetic energy. From the energy contour maps, K
+
ions experience
lower interaction energy as compared to the Cl
-
ions between adjacent ice molecules surrounding
the ice cavities. It is surmised that such lower energy barrier increases the probability of K
+
ion
motion across ice cavities as substantiated by its higher jump frequency.
The ion contributions in the transport mechanism of ice comes from the diffusion
coefficients as determined from MSD during simulation. Comparing the MSD profiles, ice
doped with ions is significantly higher than that of pure ice. In addition, the rate of ion motion as
observed by MSD is also higher along the axis of the externally applied electric field at the lower
temperature of 200K. At the higher simulation temperature of 220K, however, the influence of
the electric field is negligible. The self-diffusion coefficient of pure ice from our simulations are
in excellent agreement with that reported by experimental studies at approximately 6.30 x 10
-14

cm
2
/sec at 200K and increases with temperature to approximately 8.1 x 10
-13
cm
2
/sec at 220K.
In ice simulations containing ions, however, the presence of the electric field decreased the
diffusion coefficient by approximately 1 order of magnitude. Separating the calculation of the
diffusion coefficients of ice and ions, the motion for ions is more than an order of magnitude
higher than that of pure ice at 10
-10
cm
2
/sec. Between K
+
and Cl
-
ions, the diffusion coefficient
for K
+
ions is also substantially higher by an order of magnitude. From the diffusion

171

coefficients, the calculated conductivity from the simulations is approximately 10
-7
S/cm at 220K
and decreases by an order of magnitude at 200K. This is in excellent agreement with our
experimental results measuring the electrical properties of ice, and thus, provides further
validation on the MD simulation results on the transport mechanisms of ice. In conclusion, the
ion transport within the ice lattice can be best understood as a complex process initiated by local
melting of ice structures near the migration ion followed by frequent events of ion hopping
across adjacent ice cages. The influence of the electric field appear to be more significant at
lower temperatures as kinetic energy is more dominate at higher temperatures.

172

Chapter 5: Final Conclusions
This dissertation accomplishes the three primary goals: 1) the design and validation of an
in-situ remote sensing platform using electrochemical impedance spectroscopy technique (EIS),
2) obtain an comprehensive experimental database of the electrical properties of pure and ion
doped ice to support future fieldable instrument development, and 3) perform a computer
experiment based on atomistic simulations in an effort to provide molecular details on the
electrical properties of pure and doped ice from EIS measurements.
Chapter 2 addresses the experimental design. The probe configuration for experimental
setup using EIS is specifically tailored to grow and to perform in-situ measurements on ice. The
sense probe configuration incorporates a 32 ml Teflon cup to house the sample material with a
set of contact polarizable electrodes located at the bottom. The footprint of sense electrodes is
very portable at approximately 1 cm
2
and can be easily scalable to larger field applications.
Calibration of the measurement system indicates excellent sensitivity down to 1 mM
concentration. Using a fully automated experimental setup with controls over both the thermal
chamber and the measurement instrument, a substantial experimental database was obtained over
a wide range of temperatures (down to -65
o
C) on doped ice materials from simple electrolyte
solutions containing KCl, NaCl, MgSO
4
, KOH, and CaCl
2
.
From Chapter 3, we characterized the experimental measurements on ice and revealed
numerous electrical properties are unique to its molecular constituents and structure. From the
relative permittivity spectra, the well-known Debye relaxation was readily observed as the slow
transition between the high frequency infinite relative permittivity (ε
∞
) and the low frequency
static relative permittivity (ε
s
). The Debye relaxation has a strong time and temperature
dependency enabling distinctive detection and characterization. Unique to our electrode
configuration, resonance was also detected occurring at the onset of the Debye relaxation. The
observed resonance is attributable to the electric field orientation from the unique probe
configuration and can also be employed as a detection characteristic on the electrical properties
of ice.
By far, the conductivity measurement on ice is the most viable electrical property for ice
characterization. Ice conductivity is dependent on temperature, ion concentration, and ionic

173

species. The measured conductivity of different doped ices can exceed three orders of
magnitude higher than that of pure ice, and suggests water-like behavior with free ion motion
rather than immobile ion locked in the ice cavities. The doped ice from KCl and NaCl
electrolytes exhibits the very similar conductivity profile of pure ice due to its low ion
concentration. From visualization of these ices in polystyrene cup, we were able to observe that
the majority of ionic species are insoluble in ice and precipitate toward the center of the
container during nucleation. However, the conductivity of KCl and NaCl doped ice is still more
than an order higher than that of pure ice which is suggest a significant presence of ion
concentration into ice lattice. Ice doped with MgSO
4
exhibited nearly identical conductivity as
that of pure ice indicating virtually no ion concentration.
Modified equivalent circuit models were developed based on the Debye circuit to
effectively model both Debye relaxation and resonance behavior. From these equivalent circuit
model parameters, activation energy of pure ice was determined at 0.579 eV and is in excellent
agreement with the baseline parallel plate electrode configuration of 529 eV. With the exception
of MgSO
4
doped ice, all the activation energies of the conductivity, Debye relaxation and
resonance times for ice doped with various ions are significantly lower but reached a relatively
constant value beyond 100 mM ion concentration. We attributed this to the ion concentration in
the ice as well. Using known solubility standards at liquid conditions, we were able to estimate
this ion saturation concentration in the KCl-doped ice to be approximately 2 μM for ice.
Compared with pure ice, the estimated ion saturation concentration for all the ionic species
varies appreciably which provides another viable detection parameter in ice. The estimated
diffusion coefficient of doped ices is significantly higher than expected at 10
-5
cm
2
/sec which
indicates polycrystalline as oppose to single crystalline ice. Overall, our experimental
measurements on doped ices using EIS demonstrates an unambiguous measurement of the
electrical properties in ice which are directly attributable to ion concentration and ionic species.
From the classical understanding, our experimental results appear to support the theory that
certain ionic species are soluble in ice and can ultimately result in the increase of protonic
defects (e.g., Bjerrum and ionic types) which directly lead to the observed increases in the
transport properties. However, the true nature in which the role of ions plays on defect formation
and migration mechanisms cannot be ascertained from our experimental efforts.

174

Finally, Chapter 4 of this investigation focuses on developing a molecular interpretation
and understanding of the electrical property measurements on ice through the use of atomistic
simulation techniques. To full realize this, we performed large-scale simulations over a very
long duration of greater than 1 μ-sec in order to capture the critical structural and transport
properties in ice. Initial evaluation of ice lattices using numerous water force-field models
(TIP4P, TIP4P-Ew, TIP4P/Ice, and TIP4P/2005) revealed all water model produces the accurate
ice physical properties, including the density of about 920 Kg/m
3
and a heat capacity of about
50.0 J/K-mole. We developed large production ice lattices containing more than 6144 TIP4P
water molecules with the cubic lattice dimension of 53.9 Å x 62.2 Å x 58.7 Å. We incorporated
K
+
and Cl
-
ions and an applied electric field in our MD simulations in order to mimic actual
experimental conditions. Simulations were conducted on three unique initial ion configurations
to ensure experimental integrity and repeatability of results.
Our molecular dynamic simulations on ice revealed some critial microscopic details on
structural and transport dynamics not obtainable from macroscopic experiments such as EIS.
The RDF of K-Cl, K-K, and Cl-Cl indicated significant structural changes attributable to ion
motion throughout the simulations, but the overall ice lattice structure remained structurally
stable even under the influence of the electric field. A comparison of RDF of H-K
+
and H-Cl
-

also indicates that Cl
-
ions are more closely packed to the ice molecules which indicate Cl
-
ions
are more immobile than K
+
ions once solvated into the ice lattice. An important finding is that
ions in the ice lattice exhibits migration in the form of erratic hops across adjacent ice cavities.
Such ionic jumps were also substantiated by the trajectory plots of ion positions throughout the
simulation.
To quantify jump frequencies and timescale and offer details on the jump mechanisms,
we calculated the individual ion displacement profiles and determined the jump rates using a
technique defined by the standard deviation of the displacement running average over short
displacement period. We found the jump rate of K
+
ions is higher by a factor of approximately 5
over the jump rate of Cl
-
ions. As expected, higher frequency of jump events occur with the
increase in temperature from 200K to 220K by a factor of 3-4 due to the added thermo energy.
The timescale of such ion jump events appears to be randomly distributed over the entire 1 μ-sec
simulation time. However, the high probability of jump frequencies at the relatively short

175

timescale of less than 10000 ps suggests that ion jump events occur more frequently after the
onset of the initial jump event. From the local potential energy profile between the surrounding
ice and the ions, we found a direct correlation as the local potential energy decreases after jump
events. An energy contour analysis of the K
+
ion further indicated lower interaction energies
with the surrounding ice cages thereby causing wider migration pathways than the Cl
-
ion.
Ostensibly, the local ice structure surrounding the ions does plays a substantial role in the overall
ion transport in ice as it is determined that these local water structures are more disordered
exhibiting behaviors similar to liquid-like water. These ice structures can be also construed as
crystalline defects which seem to support the classical thought of ice transport mechanisms based
on protonic-type defects. However, this study also found that the ions themselves can contribute
significantly to the ice transport properties in such a mechanism as can be described as hopping
across ice cavities. Once initiated, such ion transport processes in the ice are highly temperature
dependent and appears to occur rapidly. From trajectory plots targeting the closer examination
of local structures near the individual ions, we were able to confirm the existence of such random
ion jump across adjacent ice cavities by both K
+
and Cl
-
ions as well as highly disordered ice
structures adjacent to the ions. The overall phenomenon is quite complex but the ion transport
processes are best be interpreted from observations of our simulations is the temporary formation
of disordered ice structures most likely due to local ice melting; these structures have lower
migration energy barriers which can facilitate frequent and rapid ion motion across adjacent ice
cavities.
To directly compare with experimental results, we also calculated the electrical properties
of doped ice using the simulation data. The diffusion coefficients were determined using the
MSD profiles. The diffusion coefficient of pure ice from simulation is in the order to 10
-13

cm
2
/sec and is in good agreement with experimental values even at 220K. By decoupling the
ions and ice molecules separately within the same simulated system, it is determined that the ice
molecules exhibit substantially lower MSD as compared to the average ion. The estimated
diffusion coefficient of the average K
+
ion is in the order of 10
-9
cm
2
/sec, nearly 1 order of
magnitude higher than that of pure ice. As expected, the diffusion coefficient is also dependent
on both temperature and the external electrical field but the influence of the electric field is
stronger at lower temperatures due difference in kinetic energy. The electrical properties are also
significantly dependent of the presence of ions. In pure ice, the dielectric constant from

176

simulations is relatively temperature independent similar to experimental results. In ice with
ions, however, both the dielectric constant and conductivity increase with increasing
temperature. Overall, the calculated electrical properties of both dielectric constant and
conductivity at the end of simulation are in excellent agreement with our experimental results at
2.39 and 10
-7
S/cm, respectively.
5.1 Suggested Future works
The experimental and atomistic simulation work represented in this study is one of the
most comprehensive investigations on physics of ice. Ice, as it exists in nature, exhibits
numerous unique physical and chemical properties contributing to its observed transport
mechanism. From our simulation work, we now have a better understanding on the ice electrical
properties. However, since experiments have shown that the transport properties of ice varies
significantly with different ionic species additional simulations comparing the transport
properties of ice doped with a variety of other ions would complement this study well. More
importantly, our unique experimental setup not only validates ice characterization using IS
technique but also provides an in-situ design platform for a portable and fieldable
instrumentation development. However, future efforts must also focus on the development of a
portable brass-board instrumentation platform suitable for field testing. Therefore, more
simulation work probing the interfacial properties of ice with other bulk materials such as soils
needs to be conducted in order to fully characterize more realistic experimental conditions with
heterogeneous mixtures.

177

Appendix A: Fundamental Equations for Conduction &
Polarization by AC Impedance

AC impedance, more formally known as Electrochemical Impedance Spectroscopy (EIS),
enables the measurement of impedance or admittance and phase angle as a function of frequency
for the material of interest. The admittance (Y) is related to the impedance (Z) by
" '
1
jY Y jB G
Z
Y      A.1
Y
’
and Y
”
are also referred to as the conductance (G) and susceptance (B) of the response. The
admittance is related to the complex conductivity and complex permittivity by
) ( '
" ' '
      j j j
Geo
Y
o o
    A.2
where Geo is the geometrical factor of the electrode configuration. Therefore, the loss tangent
(tan δ) of the dielectric loss for ac impedance measurements is
'
"
'
'
tan


 

  
o
A.3

The loss tangent is commonly used for assessing the energy dissipation due to the electric field.
For an applied alternating potential (V
ac
), the current (I
ac
) of the response is related to the
complex conductivity by
ac ac
GeoV j I ) (
" '
    A.4
where σ” = ωε
o
ε’. The current density is obtained in terms of the complex conductivity and
complex permittivity and the electrical field by
E j E j J
o o
) ' " ( ) (
" '
          A.5

178

Therefore, the conductivity (σ’) and relative permittivity (ε
r
) are related to complex admittance
by
Geo
Y
o
'
" '
     A.6

Geo
Y
o
r
 


"
"
  A.7

179

Appendix B: Derivation of the Debye Equivalent Circuit
from Basic Principles

With an application of the electric field, two principle interactions may occur in solids.
The first interaction result in a displacement current (i),
dt
dD
i  B.1
where D is the electric displacement defined by,
P E D  
0
 B.2
where E is the electric field, ε
0
is the permittivity of free space, and P is the polarization of the
dielectric material. The second interaction leads to purely real (DC) conductivity (σ),
E i   B.3
At a molecular level, the first interaction is due to the orientation of dipoles (orientational
polarization) in the electric field while the second is concerned with ionic vibrations (ionic
polarization).

180

Figure B.1: Time Dependence of Polarization with Application of Electric Field

Figure B.1 graphically depicted the time dependence of polarization P after the application of an
electric field. Simply put, the rate of P approaches P
s
is proportional to the difference between
them
) (
'
   

t P P P
s
B.4
This behavior is governed by first-order kinetics with a single relaxation time (τ) such that
P P
dt
t dP
s
 
) (
'
 B.5
With the application of a unit step voltage u
o
(t) form Figure B.1
) ( ) (
'
t P t u P P
o
 

B.6
Substituting equation B.6 into B.5 and perform Laplace transform of equation 5 (noting {u
o

(t)}=0 for t>0 by Dirac Delta Function)

P
s

P
P
∞
Time
t = 0

181

  } { ) 0 ( } { P
s
P
P P s
s
    B.7
where s is the Laplace variable and τ = 1/ω
o
. Since P(0) = P
∞
, solving for the Laplacian
dependent variable {P} in equation B.7
o o
s o
s
P
s s
P
P
 






) (
} { B.8
From equation B.8, the current density function can be determined using equation B.1
dt
dP
dt
dD
i   B.9
The Laplace Transform of equation B.9 is
) 0 ( } { } {    t P P s i B.10
Since P(t=0) = P
∞
, equation B.10 becomes

  P P s i } { } { B.11
Substituting equation B.8 into equation B.11 and simplify
o
o
s
s
P P P i



  
 
) ( } { B.12
The admittance of is simply
} {
} {
E
i
Y  B.13
where {E} is the Laplacian of the unit step voltage, {E} = 1/s. Substituting in equation B.12 in
equation B.13
o
o
s
s
s
P P sP Y



  
 
) ( B.14

182

Since polarization is independent of electric field strength, P
s
= ε
s
ε
o
and P
∞
= ε
∞
ε
o
. Noting s =
jw (complex form), equation B.14 becomes
o
o o
s o
j
j
j Y
 
 
   

  
 
) ( B.15
In terms of complex dielectric form, the equation to use is
o
j
Y

  B.16
Substituting equation B.15 in equation B.16 gives
o
o
s
j  

   

  
 
) ( B.17
By substituting in τ and simplify to typical complex format, equation B.17 becomes
) (
1 1
2 2 2 2








   
 

 
 
 
s
s
j B.18
Equation B.18 is the typical so-called Debye equation reported in literature for single dielectric
relaxation time constant in solid-state materials.

183

Figure B.2: Equivalent Circuit for Single Time Constant Solid-State Material

The equivalent circuit for the Debye equation is shown in Figure B.2. In the circuit, C1 is the
bulk capacitance, R
2
is the bulk resistance, and C2 is double layer capacitance at the electrodes.
Typically, C2>>C1 due to the minor interaction of the bulk solid material as compared to
surface electrode phenomena.
The admittance of the circuit in Figure B.2 is









 


2
2
2
1
2
2
1
2
2
2
1
2
2
2 1
2
1 1 C R
C
C j
C R
C R
Y




B.19
Starting from equation B.15, the admittance for the Debye model can also be simplified to the
form








 






2 2 2 2
2
1
) (
1
) (
 
  
  
 
    
o s
o
o s
j Y B.20
Comparing equation B.19 to B.20 provides the solutions to the circuit components of Figure B.2.
2 2
C R   B.21
o
C  


1
B.22
o s
C    ) (
2 
  B.23

C
2

C
1

R
2

184

Appendix C: Equivalent Circuit Model Derivations

The derivations of the three equivalent circuits used in this investigation are detailed in
sections C.1, C.2, and C.3.
C.1: Derivation of the Debye Equivalent Circuit

Figure C.1.1: The Debye Equivalent Circuit
Figure C.1.1 is the ideal Debye equivalent circuit model. The admittance (Y
2
) of the R
2
and C
2

circuit component in series is
2
2 2
1 1
R
C j Y
 

C.1.1
Taking the reciprocal and multiplying by the the complex conjugate, equation C.1.1 becomes
2 2
2
2
1 C R j
C j
Y



 C.1.2
Now, adding the R
o
and C
1
components to equation C.1.2 to get the total admittance for the
circuit
2 2
2
1
1 C R j
C j
C j Y




  C.1.3
Taking the complex conjugate the last term, equation C.1.3 becomes
2
2
2
2
2
2
2
2
2
2
2
2 2
2 2
1
1 1 C R
C j
C R
C R
C j Y








  C.1.4

185

Separating the real and imaginary parts, equation C.1.4 becomes







 


2 2
2
1
2
2
2
2
2
2 2
2 2
1 1 C R j
C j
C j
C R
C R
Y





C.1.6
To transform equation C.1.6 from complex admittance to the complex relative permittivity form,
we use the equation
Geo j Y
*
  C.1.7
where Geo is the cell constant and ε* the complex relative permittivity. Solving for ε* from
equation C.1.7, Equation C.1.6 becomes













 


2 2
2
1
2
2
2
2
2
2 2
2 2 *
1 1
1
C R j
C j
C j
C R
C R
G j
eo






 C.1.8
The circuit components of R
o
, C
1
, and C
2
can be re-defined as follows


  
 

 



2 2
2
1
) (
C R
G C
G C
eo s o
eo o
C.1.9
Substituting the components in C.1.9 into equation C.1.8 and simplify, the complex relative
permittivity equation becomes
2 2 2 2
*
1
) (
1
) (
 
  

 
 








 
 

s s
o
j C.1.10
Finally, Equation C.1.10 can be further simplified to the following generalized form

 
  


j
j
s
o


   


1
"
'
*
C.1.11

186

C.2: Derivation of the Debye Equivalent Circuit with Leakage.

Figure C.2.1: The Debye Equivalent Circuit with Leakage

Figure C.2.1 is the Debye equivalent circuit model with leakage. The admittance (Y
2
) of the R
2

and C
2
circuit component in series is
2
2 2
1 1
R
C j Y
 

C.2.1
Taking the reciprocal and multiplying by the the complex conjugate, equation C.2.1 becomes
2 2
2
2
1 C R j
C j
Y



 C.2.2
Now, adding the R
o
and C
1
components to equation C.2.2 to get the total admittance for the
circuit
2 2
2
1
1
1
C R j
C j
C j
R
Y
o




   C.2.3
Taking the complex conjugate the last term, equation C.2.3 becomes
2
2
2
2
2
2
2
2
2
2
2
2 2
2 2
1
1 1
1
C R
C j
C R
C R
C j
R
Y
o








   C.2.4
Separating the real and imaginary parts, equation C.2.4 becomes

187








 

 
2 2
2
1
2
2
2
2
2
2 2
2 2
1 1
1
C R j
C j
C j
C R
C R
R
Y
o





C.2.6

To transform equation C.2.6 from complex admittance to the complex relative permittivity form,
we use the equation
Geo j Y
*
  C.2.7
where Geo is the cell constant and ε* the complex relative permittivity. Solving for ε* from
equation C.2.7, Equation C.2.6 becomes













 

 
2 2
2
1
2
2
2
2
2
2 2
2 2 *
1 1
1 1
C R j
C j
C j
C R
C R
R G j
o eo






 C.2.8
The circuit components of R
o
, C
1
, and C
2
can be re-defined as follows


  
 


 




2 2
2
1
) (
1
C R
G C
G C
G
R
eo s o
eo o
DC eo
o
C.2.9
Substituting the components in C.2.9 into equation C.2.8 and simplify, the complex relative
permittivity equation becomes
o
DC s s
o
j j


 
  

 
 









 
 

2 2 2 2
*
1
) (
1
) (
C.2.10
Finally, Equation C.2.10 can be further simplified to the following generalized form
o
DC s
o
j
j
j



 
  





   


1
"
'
*
C.2.11

188

C.3: Derivation of the Debye Equivalent Circuit with Leakage and Resonance.

Figure C.3.1: The Debye Equivalent Circuit with Leakage and Resonance
Figure C.3.1 is the Debye equivalent circuit model with leakage and resonance. The admittance
(Y
2
) of the R
2
, L
2
, and C
2
circuit component in series is
2 2
2 2
1 1
R L j
C j Y
   

C.3.1
Taking the reciprocal and multiplying by the the complex conjugate, equation C.3.1 becomes
2 2 2 2
2
2
2
1 C R j C L
C j
Y
 

 
 C.3.2
Next, the natural resonance frequency (ω
o
) is defined as
5 . 0
2 2
) (
1
C L
o
  C.3.3
Substitute equation C.3.3 into C.3.2 and simplify to complex form, the admittance (Y
2
) becomes
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2 2
2
2
1
1
1 C R
C
j
C R
C R
Y
o
o
o








































 C.3.4

189

Now, adding the R
o
and C
1
components to equation C.3.4 to get the total admittance for the
circuit
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2 2
2
1
1
1
1
1
C R
C
j
C R
C R
C j
R
Y
o
o
o
o









































   C.3.5
Separating the real and imaginary parts, the final form of the equivalent circuit is

































 










 
2
2
2
2
2
2
2
2
2
2
2
1
2
2
2
2
2
2
2
2
2
2 2
2
1
1
1
1
C R
C
C j
C R
C R
R
Y
o
o
o
o










C.3.6
To transform equation C.3.6 from complex admittance to the complex relative permittivity form,
we use the equation
Geo j Y
*
  C.3.7
where Geo is the cell constant and ε* the complex relative permittivity. Solving for ε* from
equation C.3.7, Equation C.3.6 becomes















































 










 
2
2
2
2
2
2
2
2
2
2
2
1
2
2
2
2
2
2
2
2
2
2 2
2
*
1
1
1
1 1
C R
C
C j
C R
C R
R G j
o
o
o
o eo











 C.3.8
The circuit components of R
o
, C
1
, and C
2
can be re-defined as follows

190


  
 


 




2 2
2
1
) (
1
C R
G C
G C
G
R
eo s o
eo o
DC eo
o
C.3.9
Substituting the components in C.3.9 into equation C.3.8 and simplify, the complex relative
permittivity equation becomes
   
 
 
 









 


 
 
 
 

0
2 2
2
2
0
2 2 2
2
2
0
2
2
0
2
*
/ 1 / 1
/ 1


   
  
   
   
 
DC
D
s
D
s
j C.3.10
Finally, Equation C.3.10 can be further simplified to the following generalized form
o
DC
o
s
o
j
j
j





 
  



 

   


2
2
'
*
1
" C.3.8

191

Symbols
a
A+
Chemical Activity of species A
+
.
a
B+
Chemical Activity of species B
+
.
b Fitting parameter for Debye-Huckel extended law.
c
p
Heat capacity, J/K-mole
e Unit of charge, 1.602 x 10
-19
C.
g(r) Radial distribution function.
g(O-O) Radial distribution function between the Oxygen atoms.
g(K-Cl) Radial distribution function between the K
+
and Cl
-
atoms.
g(K-K) Radial distribution function between the K
+
atoms.
g(Cl-Cl) Radial distribution function between the Cl
-
atoms.
g(H-K) Radial distribution function between the H and K
+
atoms.
g(O-K) Radial distribution function between the O and K
+
atoms.
g(H-Cl) Radial distribution function between the H and Cl
-
atoms.
g(O-Cl) Radial distribution function between the O and Cl
-
atoms.
j Complex number, 1  , or current density.
m
i
Mass of solute.
mM Milli-Mole, 10
-3
Moles.
n Number of ionic species.
k Coefficient describing system of harmonic oscillators.
k
b
Boltzman constant, 1.38 x 10
23
J/K.
p Exponential fitting factor for constant phase element.
p Dipole moment vector, units of Debye.
r
i
Stokes Radius, measures the radius of interaction between ionic species.

192

r
o
LJ characteristic distance, Ǻ.
r(0) Initial position.
s Solubility (grams of solute/liter of solvent).
t time.
i
v

Velocity vector of species i.
z
i
Charge number of ion i.
z Displacement distance.
A Coefficient describing the vibrational state of a crystal.
B Susceptance, S.
C Concentration, Moles/liter.
C
d
Double-layer capacitance, C/V.
C
i
Concentration of species i, mole/liter.
C
o
Standard concentration, 1 mole/liter.
Di Diffusion coefficient of species i, cm
2
/sec.
Dx, Dy, Dz Diffusion coefficient, cm
2
/sec, along the X,Y,Z-axis.
E
a
Activation Energy, eV.
E
s
Electric field strength, N/C.
E, E
o
Unit of Energy or Electric Field, V/m.
E
LD
Activation energy of formation for L,D-Bjerrum defects.
F Faraday’s constant, 96,400 C/mole.
F
i
Force acting of species i.
G Conductance, S.
G
eo
Geometric Factor for electrode configuration, cm
-1
.
H Unit of inductance, Henry.
I Ionic strength, moles/kg.

193

J Current density, A/cm
2
.
K Equilibrium constant.
K
sp
Solubility product.
P Polarization.
R Gas constant, 3.314 J/K-Mole.
M Dipole moments, C- Ǻ.
N
A
Avagadro’s constant, 6.022 x 10
23
atoms/mole.
N Total number of particles in a system.
N
LD
L,D-Bjerrum defect concentration at equilibrium.
T Absolute temperature, K.
U Internal Energy.
V Volume of simulation box, Ǻ
3
.
Y Admittance, S.
Z Impedance, Ω.
α Degree of dissociation (unitless), OR polarizability, Ǻ
3
.
α
p
Thermal expansion coefficient, ppm/
o
K.
μM Micro-moles, 10
-6
M.
μ-sec Micro-seconds, 10
-6
seconds.
μ
i
Mobility of species i, cm
2
/sec-V.
ε
o
Permittivity of free space, 8.85 x 10
-12
C/N-m
2

ε
r
Dielectric constant, 78.3 for water at 25
o
C.
ε(o) Static dielectric constant.
ε
∞
High (“infinite”) frequency dielectric constant.
ε
s
Low (“static”) frequency dielectric constant.
ε’ Real part of complex dielectric constant.

194

ε” Imaginary part of complex dielectric constant.
ε
ij
LJ interaction parameter for depth-well between atoms i and j, kcal/mole.
ρ
s
Solvent density, 997.05 kg/m
3
.
κ
-1
Debye screening length, Ǻ.
κ
T
Isothermal compressibility, area/force.
ν
d
Drift velocity, m/s.
σ Conductivity, S/cm.
σ’ Real part of complex conductivity, S/cm.
σ” Imaginary part of complex conductivity, S/cm.
σ
ij
LJ interaction parameter for atomic radius between atom i and j, Ǻ.
τ, τ
r
, τ
D
Relaxation time constants.
χ
s
Susceptibility.
ω Angular frequency, 2πf.
ω
o
Resonant frequency.
ω
n
Natural resonant frequency.
χ
e
Electrical susceptibility.
γ
A-
Activity coefficient of species A
-
.
γ
B-
Activity coefficient of species B
-
.
Ω Unit of resistance, Ohms.
Λ Molar conductivity, S/cm-mole.
Λ
o
Molar conductivity at infinite dilution, S/cm-mole.
 Del operator.
<…> The average of a value.
{…} Laplacian variable.

195

Bibliography

Allen, M.P. and Tildesley, D. J., Computer Simulations of Liquids, Clarendon Press, Oxford,
1987.
Anderson, J.C., Dielectrics, Spottiswood, Ballantyne & Co Ltd, Great Britain, 1964.
Bard, A.J. and Faulkner, L.R., Electrochemical Methods: Fundamentals and Applications, 2
nd

Ed., John Wily & Sons, NY, 2001.
Barsoukov, E. and Macdonald, J. R., Impedance Spectroscopy: Theory, Experiment, and
Applications, 2
nd
Edition, John Wiley & Sons, Inc., New Jersey, 2005.
Carr, M.H., Water on Mars, Oxford University Press, NY, 1996.
Debye P., Polar Molecules, Reprinted 1945, New York: Dover, 1929.
Fenkel, D. and Smit, B., Understanding Molecular Simulations: From Algorithms to
Applications, Academic Press, London, 2002.
Fletcher, N.H., The Chemical Physics of Ice, The University Press, Cambridge, 1970.
Fröhlich H., Theory of Dielectrics, Oxford University Press, 1949.
Jakosky, B., The Search for Life on Other Planets, Cambridge: University Press, Cambridge,
1998.
Moulson, A.J. and Herbert, J.M., Electroceramics: Materials, Properties, Applications,
Chapman & Hall, London, 1992.
O’M. Bockris, J., Reddy, A.K.N., and Gamboa-Aldeco, M., Modern Electrochemistry 2A:
Fundamentals of Electrodics, 2
nd
Ed., Kluwer Academics/Plenum Publishers, New York, 2000.
Onsager, L. and Dupuis, M., Electrolytes, Pergamum Press, Inc., New York, 1962.
Pang, T., An Introduction to Computational Physics, 2
nd
Ed., Cambridge University Press,
Cambridge, 2006.
Petrenko, V.F. and Whitworth, R.W., Physics of Ice, Oxford University Press, 1999.
Rowley, R.L., Statistical Mechanics for Thermophysical Properties Calculations, PTR Prentice
Hall, New Jersey, 1994.
Van Valkenburg, M.E., Network Analysis, 3
rd
Ed., Prentice-Hall, Inc., NJ, 1974.
von Hippel, A., Dielectrics and Waves, J. Wiley, New York, 1954.

196

References

1
E. J. Gaidos, K. H. Nealson, J. L. Kirschivnk, Science, 284 1631 (1999).
2
C.D. Parkinson, M.C. Liang, Y. L. Yung, J. L. Kirschivnk, Astrobiology, 38, 355 (2008).
3
D. Schulze-Makuch, J. M. Dohm, A. G. Fairen, V. R. Baker, W. Fink, & R. G. Strom,
Astrobiology, 5, 778 (2005).
4
S. W. Squyres et al., Science, 306, 1709 (2004).
5
O. Botta, J. L. Bada, J. Gomex-Elvira, E. Javaux, F. Selsis, R. Summons, Space Science
Review. 135, 371 ((2008).
6
B. M. Jakosky, K. H. Nealson, C. Bakermans, R. E. Ley, & M. T. Mellon, Astrobiology, 3, 343
(2003).
7
M. Cabane, et al., Advances in Space Research, 33, 2240 (2004).
8
G. R. Olhoeft, Geophysics, 50, 2492 (1985).
9
M. G. Buehler, K. B. Chin, S. Seshadri, D. Keymeulen, R. C. Anderson, T. A. McCann, IEEE
Aerospace Conference Proceedings 2007, Piscataway, NJ, doi 10.11109/AERO.2007.352776.
10
S. Seshadri, K. B. Chin, M.G. Buehler and R.A. Anderson, Astrobiology, 8, no. 4, 781 (2008).
11
R.P. Auty and R. H. Cole, Journal of Chemical Physics, 20, 1309 (1952).
12
G.P. Johari and E. Whalley, Journal of Chemical Physics, 75, 1333 (1981).
13
O. Wörz and R. H. Cole, Journal of Chemical Physics, 51, 1546 (1969).
14
R.E. Grimm, Journal of Geophysical Research, 107, doi 10.1029/2001JE001616, (2002).
15
J. Fleig and J. Maier, Electrochemica Acta, 41, 1003 (1996).
16
Y. Feldman, R. Nigmatullin, E. Polygalov, J. Texter, Physical Review E, 58, 7561 (1998).
17
R. Hassan and E.S. Campbell, Journal of Chemical Physics, 97, 4326-4335 (1992).
18
P.L.M. Plummer, Physics and the chemistry of ice, Hokkaido University Press, Sapporo, 54-61
(1992).
19
R. Pedeszwa and V. Buch, Physical Review Letters, 83, no. 22, 4570 (1999).
20
N. Grishina and V. Buch, Journal of Chemical Physics, 120, no. 11, 5217 (2004).

197

21
S. Kawada, H. Iisada, E. Kitamure, R. Tutiya, and H. Abe, Physics and the chemistry of ice,
Hokkaido University Press, Sapporo, 20-26 (1992).
22
S. Kawada, Journal of the Physical Society of Japan, 58, no. 1, 54 (1989).
23
S. Kawada, R.G. Jin, and M. Abo, Journal of physical Chemistry, B101, 6223-6225 (1997).
24
J.D. Bernal and R.H. Fowler, Journal of Chemical Physics, 1, 515 – 548 (1933).
25
N. Bjerrum, Nature, 385 (1952).
26
L. Pauling, Journal of the American Chemical Society, 57, 2680 (1935).
27
J.G. Kirkwood, Journal of Chemical Physics, 4, 592 (1936).
28
J.G. Powles, Journal of Chemical Physics, 20, 1302 (1952).
29
M. Neumann, Molecular Physics, 4, 841-858 (1983).
30
C. P. Smyth and C. S. Hitchcock, Journal of American Chemical Society, 54, 4631 (1932).
31
E.J. Murphy, Transactions of Electrochemical Society, 65, 309 (1934).
32
R.P. Auty and R.H. Cole, Journal of Chemical Physics, 20, 1309 (1952).
33
H. Granicher, Proceeding of Royal Society, A247, 453 (1958).
34
C. Haas, Physics Letters, 3, 126 (1962).
35
L. Onsager and L. K. Runnels, Journal of Chemical Physics, 50, 1089 (1969).
36
B. Geil, T.M Kirschgen, F. Fujara, Physical Review B, 72, 014304 (2005).
37
J.F. Nagle, Chemical Physics, 43, 317 – 328 (1979).
38
A. von Hippel, The Journal of Chemical Physics, 54, 145 (1971).
39
C. Jaccard, Helvetica Physica Acta, 32, 89 (1959).
40
A. von Hippel, D.B. Noll, and W.B. Westphal, The Journal of Chemical Physics, 54, no. 1,
134 (1971).
41
M.A. Maidique, A. von Hippel, and W.B. Westphal, The Journal of Chemical Physics, 54, no.
1, 150 (1971).
42
S.R. Gough and D.W. Davidson, Journal of Chemical Physics, 52, no. 10, 5442 (1970).
43
G.C. Camplin, J.W. Glen, and J.G. Paren, Journal of Glaciology, 21, No. 18, 123 (1978).

198

44
G.P. Johari and E. Whalley, Journal of Chemical Physics, 115, 3274 (2001).
45
V.F. Petrenko and I.A. Ryshkin, Journal of Physical Chemistry, B101, 6285 (1997).
46
G.P Johari and S.J. Jones, Journal of Glaciology, 21, 259-276 (1978).
47
S. Kawada, Journal of the Physical Society of Japan, 44, 1881-1886 (1978).
48
A. von Hippel, Journal of Electrochemical Society, 199, no. 2, 45C (1972).
49
E.J. Workman and S.E. Reynolds, Physical Review, 78, no. 3, 254 (1950).
50
A.V. Zaretskii, V.F. Petrenko, A.V. Trukhanov, E.A. Aziev, and M.P. Tonkonogov, Journal de
Physique, 48, 87-91 (1987).
51
K. Goto, T. Hondoh, and A. Higashi, Japanese Journal of Applied Physics, 25, no. 3, 351-357
(1886).
52
A. K. Soper, Chemical Physics, 258, 121-137 (2000).
53
J.L.F. Abascal, R.G. Fernández, L.G. MacDowell, E. Sanz, C. Vega, Journal of Molecular
Liquids, 136, 214–220 (2007)
54
C. Vega, C. McBride, E. Sanz, and J.L.F. Abascal, Physical Chemistry Chemical Physics ,7 ,
1450–1456 (2005)
55
E.J. Smith and A.D.J. Haymet, Molecular Simulation, 30, No. 13-15, 827–830 (2004).
56
P.J. Feibelman, Physical Review B, 75, 214113, (2007).
57
D.J. Adams, J. Phys. C: Solid State Phys., 17, 4063-4070 (1984).
58
H. Itoh, K. Kawawmura, T. Hondoh, S. Mae, Physica B, 219-220, 469-472 (1996).
59
T. Ikeda-Fukazawa, S. Horikawa, and T. Hondoh, The Journal of Chemical Physics, 117, no.
8, 3886-3896 (2002).
60
T. Ikeda-Fukazawa, K. Kawamura, and T. Hondoh, Molecular Simulation, 30, 15 November–
30 December, 973–979 (2004).
61
S.G. Kalko, G. Sese, and J.A. Padro, Journal of Chemical Physics, 104, No. 23, 9578-9585
(1996).
62
A.P. Lyubartsev, O.K. Forrisdahl, and A. Laaksonen, J. Chem. Phys. 108 (1), 227-233 (1998).
63
T. Bryk, A.D.J. Haymet, Journal of Molecular Liquids, 112, 47–50 (2004).

199

64
E.J. Smith, T. Byrk, and A.D.J. Haymet, The Journal of Chemical Physics, 123, 034706
(2005).
65
M.A. Carignano, P.B. Shepson, I. Szleifer, Molecular Physics, 103, 2957-2967 (2005).
66
J.S. Kim and Arun Yethiraj, The Journal of Chemical Physics, 129, 124504, (2008).
67
X.F. Pang, European Physics Journal B, 49, 5–23 (2006).
68
W. Evans, J. Fish, and P. Keblinski, The Journal of Chemical Physics, 126, 154504 (2007).
69
J.Y. Yan and G.N. Patey, J. Phys. Chem. A, 116, 7057−7064 (2012).
70
D.H. Jung, Jung Hwan Yang, Mu Shik Jhon, Chemical Physics, 244, 331–337 (1999).
71
S. Wei, C. Zhong, and H. Su-Yi, Molecular Simulation, 31, No. 8,, 555–559 (2005).
72
K. Kadota, A. Shimosaka, Y. Shirakawa and J. Hidaka, Journal of Nanoparticle Research, 9,
377–387(2007).
73
A. Vegiri and S. V. Schevkunov, The Journal of Chemical Physics, 115, no. 9, 4175-4185
(2001).
74
A. Vegiri and S. V. Schevkunov, The Journal of Chemical Physics, 116, no. 20, 8786-8798
(2002).
75
S.W. Rick and A.D.J. Haymet, Journal of Chemical Physics, 118, no.20, 9291-9296 (2003).
76
Y. W. Tang, Q. Zhang, K.Y. Chan, Chemical Physics Letters 385, 202–207 (2004).
77
F. Noritake, K. Kawamura, T. Yoshino, E. Takahashi, Journal of Non-Crystalline Solids, 358,
3108-3118 (2012).
78
K. Kadota, A. Shimosaka, Y. Shirakawa and J. Hidaka, Journal of Nanoparticle Research, 9,
377–387(2007).
79
A. Rahman, Physics Review, 136, A405-A-411 (1964).
80
L. Verlet, Physical Review, 165, 201-214 (1968).
81
D. Levesque and L. Verlet, Physical Review A, 2, 2514-2528 (1970).
82
C. Calemen and David van der Spoel, Angew. Chem., 120, 1439-1442 (2008).
83
M. Warren and J. Rottler, EPL, 88, 58005 (2009).
84
W. L. Jorgensen, J. Chandrasekhar, and J. D. Madura, J. Chemical Physics, 79, 926-935
(1983).

200

85
H. J. C. Berendsen, J.R. Grigera, and T. P. Straatsma, The Journal of Physical Chemistry, 91,
no. 24, 6269-6271, (1987).
86
J. Aqvist, J. Physical Chemistry, 94, 8021-8024 (1990).
87
D.E. Brown and S.M. George, Journal of Physical Chemistry, 100, 15460-15469 (1996).
88
R.O Ramseier, Journal of Applied Physics, 38, 2553-2556 (1967).

Abstract (if available)

Abstract As mankind approaches the cusp of deep space exploration in our quest to seek out new or similar forms of life, ice has experienced a resurgence of new research interests due to its omnipresence beyond Earth. Therefore, the search for life beyond our planet hinges on the understanding of planetary ice and its interactions with bio‐precursors that could preserve indications of pre‐existent life or currently capable of sustaining ecosystems suitable for life by providing nutrients in the form of chemical energy. ❧ The dissertation manuscript investigates the electrical properties of ice using electrochemical impedance spectroscopy (EIS) and atomistic simulations in an effort to develop an unambiguous instrumentation platform for in‐situ planetary exploration. The experimental portion of the work emphasized the electrical properties behavior of ices doped with ion in order to characterize unique the signatures from detectable from EIS measurements on different ionic species. To further delve into the transport properties of pure and ion doped ices at the molecular level, molecular dynamics simulations were also performed on KCl‐doped ice under various temperature conditions and over a long simulation time of greater than 1 μ-sec. The systems of atomistic simulations on ice included the addition of the electric field to mimic actual experimental conditions. ❧ The experiment process also included the development of a unique electrode configuration suitable for long‐life field applications. Measurement results on pure ice using this novel electrode configuration are in excellent agreement with measurements from baseline parallel plate electrodes. Resonance is also detected with the new configuration at the onset of the Debye relaxation which is unique to our measurement results. A database of ice measurements resulted in a substantial collection of ice electrical properties parameters suitable for ice characterization based on ionic species and concentration. ❧ From simulations on pure and KCl‐doped ices, structural and transport dynamics was characterized using RDF, MSD, and displacement profiles. In addition, the physical properties (e.g., density and heat capacity) of ice calculated from simulation are in agreement with actual ice indicating good simulation behavior and provide further validation of results. Comparison of the electrical properties from the simulations and the experiments indicated excellent agreement as well. Though many new finding were made regarding the dynamic behavior of ice, the core discovery from simulation is the motions of ions across adjacent ice cavities due to erratic jumps which enhance the observed transport properties. Such jump events were also substantiated by trajectory plots over the entire simulation time.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Molecular dynamics studies of protein aggregation in unbounded and confined media

PDF

Molecular dynamics simulation study of initial protein unfolding induced by the photo-responsive surfactants, azoTAB

PDF

Experimental study and atomic simulation of protein adsorption

PDF

Atomistic simulation of nanoporous layered double hydroxide materials and their properties

PDF

Molecular dynamics study of energetic materials

PDF

A molecular dynamics study of interactions between the enzyme lysozyme and the photo-responsive surfactant azobenzene trimethylammonium bromide (azoTAB)

PDF

Dynamics of water in nanotubes: liquid below freezing point and ice-like near boiling point

PDF

Experimental and atomistic simulation studies of water sorption in polyaniline

PDF

Growth, characterization of gallium arsenide based nanowires and application in photovoltaic cells

PDF

Kinetic studies of collisionless mesothermal plasma flow dynamics

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Numerical and experimental investigations of ionic electrospray thruster plume

PDF

The cathode plasma simulation

PDF

Particle-in-cell simulations of kinetic-scale turbulence in the solar wind

PDF

Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

PDF

Interactions of planetary surfaces with space environments and their effects on volatile formation and transport: atomic scale simulations

PDF

Investigation of the electrode-tissue interface of retinal prostheses

PDF

Grid-based Vlasov method for kinetic plasma simulations

PDF

Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations

PDF

Electrical signals to assess cardiovascular phenomenon

Asset Metadata

Creator Chin, Keith (author)

Core Title Ice: an impedance spectroscopy and atomistic simulation study

School Andrew and Erna Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Chemical Engineering

Publication Date 07/01/2014

Defense Date 05/22/2014

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag atomistic simulation,Ice,impedance spectroscopy,OAI-PMH Harvest,protonic defects

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Shing, Katherine (committee chair), Gupta, Malancha (committee member), Nakano, Aiichiro (committee member)

Creator Email keith.b.chin@jpl.nasa.gov,keithchi@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-427615

Unique identifier UC11286569

Identifier etd-ChinKeith-2596.pdf (filename),usctheses-c3-427615 (legacy record id)

Legacy Identifier etd-ChinKeith-2596.pdf

Dmrecord 427615

Document Type Dissertation

Format application/pdf (imt)

Rights Chin, Keith

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