This model calculates changes in air mass properties following their outflow (with or without cooling and/or windpeed and/or stable O isotope anomalies) over the sea, where evaporation adds vapour (with its own stable O composition that derives from the sea surface properties). The evaporation lowers the temperature (SST) the surface-water mixed layer (depth of which is set with the property "mld" - for temperature change we also account for insolation "Ins"), and the freshwater removal by evaporation also changes its salinity (SSS). Evaporation is calculated in daily increments, and the travel time over water is determined by means of the travel distance ("Le" for Israel and "Le/2" for Crete). We calculate all air- and water-mass properties and changes over a unit area of 1 m2. For air properties at the end-point of the trajectory to Crete, we use spacC=1, and for Haifa (Israel) we use tuning factor spacH, as explained in the main text. These tuning factors are set once, fot ensure that - in control scenarios - the Crete and Haifa relative humidities and precipitation values match with observations. Then, we keep everything constant with respect to "spac". The "spac" factor in effect stands for the lateral (or vertical) spreading of air masses as they move away from the North Aegean. But more importantly, it allows tuning of the model to modern observations with one single parameter.
Total pathway over eMed from N Aegean to Israel is roughly 1400 km
For N Aegean to Crete we use half that.
Wind speed (double peak annually in mid winter and in mid summer) after Santorini: data on https://weatherspark.com/averages/32207/Santorini-Thira-South-Aegean-Greece
mean incoming air Temp
Based on approx Plovdiv, Bulgaria, mean monthly T range from 1 to 23 deg C, with annual mean of 12.
http://www.plovdiv.climatemps.com
seasonal inc air T range
Relative humidity
mean incoming relative humidity
Plovdiv, Bulgaria annual range is 0.45-0.75
We
seasonal inc rel humidity range
(Alexandropolis, Greece range is 0.55-0.75 acc to climatemps
set seasonal cycle for incoming rel. humidity:
Insolation
mean incoming insolation at top of atmosphere (TOA)
seasonal inc insolation range
set a (=1-
For Athens, received irradiation goes from about 50 to 225 W/m^2, acc. to http://www.mdpi.com/2071-1050/6/8/5354/htm
mean Sea Surface Temp
T range from Rohling et al (2002 Clim Dyn), as 13 to 23 deg C.
seasonal SST range
We set D to cycle between about 1000 and 1400 m (see main text).
Santorini T range is 12-27 degrees.
http://www.climatedata.eu/climate.php?loc=grxx0030&lang=en
Heraklion T range is 15-29 degrees.
seasonal inc air T range
Haifa range is 14-28 deg C according to climatemps
seasonal inc air T range
We now import a file in which the same sinusoidal background signals (as outlined above) are perturbed with 1, 3, or 5 ten-day events in winter. The file starts with age, and then lists TaO (for events between -0 and -20 deg C), V (with events between 0 and 20 m/s), TaC0 for Crete (with events between -0 and -20 deg C), TaH0 for Haifa
Read out the various parameter files from "M"
References in addition to those in the main text:
Garrett, C., Outerbridge, R., Thompson, K., 1993. Inter-annual variability in Mediterranean heat budget and buoyancy fluxes. J. Climatol. 6, 900-910.
Mparsakis, S. I., 1999, Synoptic and Thermodynamic Analysis of Snowfall in the Region of Thessaloniki, Master Thesis, University of Thessaloniki (in Greek).
Wells, N.C., 1995. Surface heat fluxes in the western equatorial Pacific Ocean. Ann. Geophys., 13, 1047-1053.
Wells, N.C., 1986. The Atmosphere and Ocean: a physical introduction, 347pp. Taylor and Francis, Bristol, Pa.
Then, mean wind speeds over distances to Crete and Haifa are:
set selector for "height" so that the value "maxhei" is returned for the cold events. Because of rounding, the two critical terms are not exactly equal, so need to set a threshold of one tenthousandth (maxheiC and maxheiH are set below on the settings-page, page 9).
Start with those over the Aegean sector (Northern Greece to Crete)
Notes to p.4:
Values from page 41 (citing Korres et al., 2000) of https://books.google.com.au/books?id=w0bVBgAAQBAJ&pg=PA34&lpg=PA34&dq=aegean+and+albedo&source=bl&ots=QAGgvMZ-FQ&sig=65GJ-8H1L-tAa0e9RK4AF0IAur0&hl=en&sa=X&ei=XBEnVe7WOM6F8gWN-YDQDA&ved=0CE4Q6AEwCA#v=onepage&q=aegean%20and%20albedo&f=false
Use lapse rate to calculate air T change due to uplift over Crete
over Haifa (Israel)
and to the mean height of air column D:
Point in proof:
Gat et al. (2005) works through an adjusted "r" for the SST. This method gives the exact same evaporation as my method.
Gat et al. argue to calculate vapour pressure at incoming T, and then calculate how much moisture is imported (M0). Then RECALCULATE relative humidity for the SST, because that's the T at which evaporation would take place. That rel hum will be MUCH lower than the imported rel. hum. (but of course one then needs to use qs rather than qa in the evap calculation). If the model setup is correct, this will give identical answers:
Calculate saturation equation at sea surface temperature
Calculate saturation equation at air temperature (mean of D) during evaporation
The solutions are identical:
Watt (= J/s) per m^2
Where the latter equation (LH) gives the Latent Heat
4a
specific heat capacity of water:
See http://www.jgsee.kmutt.ac.th/exell/JEE661/JEE661Lecture2.html
no. of seconds in a day:
LH*m gives the daily heat loss in Joules/m^2
We are looking at a loss per m^2 from a surface layer of "mld" thickness, i.e., from mld m^3 of water, which is about 1.040*mld*10^3 kg
where "mld" is the mixed layer depth set in the starting conditions
Daily change in surface water T due to evaporative heat loss:
First we look at the annual T loss due to evaporation
in the 'control' years (i.e., 1-event years):
In measurements, seasonal SST range is about 10 degrees (13 to 23 deg C -see above) - and our Evaporation-caused T loss must be approx. similar
Next, the total annual T change (after insolation warming is added) is:
So the annual T loss from an "mld" m layer of water is given, just from evaporative cooling. The observed Aegean range of 10 deg C is closely found when mld = 100. Clearly, evaporative heat loss explains most of the seasonal cycle of T in the surface waters of the Aegean. SST cycles back up due to irradiation. So we can satisfactorily view the seasonal cycle this way, and only calculate the SST change anomalies associated with perturbations to the regular seasonal evaporation cycle, due to imposed polar air outbreaks.
By bringing in realistic insolation values, the annual net change is small - just as it should be (seeing that we are not accounting for conductive heat loss, or heat advection). Clearly, the vast bulk of the energy for evaporation comes from the insolation (as it should be).
In the "event" years, we find:
So we find annual SST change anomalies of the order of:
degrees C
psu
Bear in mind that, in SST, the events in proxy data seem to be of the order of 1 to 2 deg C, so we would want to see approx agreement here.
We are adding flux E (in daily proportions, so E/365) into air column of depth D, temperature Ta0, and relative humidity r. Consequently, r
Travel time over water:
Crete =
Haifa =
The 1/
in m/day
addition to the air, by weight (per unit area)
in kg of freshwater per m^2
r0*qa is expressed in (10^-3) kg water per kg of air
So in a colum of thickness D, there is a mass of moisture per m^2 as determined D*r0*qa/Dmean(weight of air per m^3):
Density of air at sea level =
Work out mean pressure and density for D
where H=7.4 km
http://acmg.seas.harvard.edu/people/faculty/djj/book/bookchap2.html
Calculate saturation state for Crete/Santorini Temps
Calculate saturation state for Haifa Temps
These are needed to calculate saturation pressures at mean Temp of column D
Because we assume specific humidity is constant through D, initial moisture contained in D is:
and Max moisture contained at the end-points (C and H) is:
Now add the daily evap over Tt timesteps, and then can calculate the saturation state again. We will assume that up to Tt*dayEadw is added, OR less until r endvalue (rend) =1
This is at the surface, at Ta0. At higher levels T is lower and so rend (as r) is higher, even if air and vapour are well mixed, because of lapse rate, and because we assume constant specific humidity in the column of thickness D!
max moisture potential for Crete air
Actual moisture at end of passage
Note that traveltime is in DAYS!
Relative humidity of the air (at Temp of meanD) at Crete:
For Crete:
Fraction of remaining atm. vapour remaining:
MmaxCout is the max moisture load at T after dry adiabatic descent
evaporation:
in m/day
in kg of freshwater per m^2
max moisture potential for Haifa air
Actual moisture at end of passage:
rel. hum of air arriving at surface Haifa (rendH):
For "spacH", see parameter settings page (9)
continued on page 10
mixed layer depth
Vinc up to 10 would even be reasonable as Bora winds can be consistent at 17 m/s http://www.istrianet.org/istria/geosciences/meteorology/winds/bora-adriatic.htm
Cooling may be set up to 9-10 deg C for an event as the wind entering is feezing (Casford et al., 2003). Below we argue why 6 seems the most representaitve, given the isotope results versus isotope data.
This uses Lawrence's formula (T in deg. C) for hLCL, the lifting condenstation level height (cloud base)
M. G. Lawrence, "The relationship between relative humidity and the dew point temperature in moist air: A simple conversion and applications", Bull. Am. Meteorol. Soc., 86, 225-233, 2005
sea mixed layer depth
Lapse rate over Israel typically between 6.5 and 7 (see main text)
During events, the cloud base is seen to drop to sea level, suggesting that widespread fog/icy frost conditions mat prevail.
SPATIAL expansion of column D, via an 'apparent' height calculation with a multiplier
This means we are adding 1/spac times the moisture to our column D through evaporation
Need to tune this manually (once!)
See also https://en.wikipedia.org/wiki/Equilibrium_level
We use spacH as a tuning factor for the control years, and then keep it constant while doing scenarios that look at relative changes in the event years. Note that wind speeds are proportionally changed (in a linear manner) to ensure conservation of mass.
Max height of convection is given by Equilibrium Level (or Limit of Convection). Usually is ~ the tropopause; Here we just use the observed cloud top of 5 or 6.5 km.
9a
Initial surface sea-water
should be close to 1.5 permil, after https://data.giss.nasa.gov/o18data/
continued from p.8
Fraction of remaining atm. vapour:
The difference between Mend and M1 both for C and for H is the amount of precipitation that has to follow uplift, in kg/m^2. Thus we can calculate P in mm for each region by considering that weight in mm of P over a sq m.
1 kg/m^2 = 1 mm/m^2
Precip amounts (in mm per m^2) then are:
Annual rainfall for Heraklion (PC) is 500 mm (see link above with T data) (of which 70 in summer)
Annual rainfall
Crete
Annual P
control:
Haifa
Annual P
control:
The blue highlighted numbers are tuned with hme and spacH - we set things up so we get a reasonable match for the 1-anomaly case
To determine E-P in Crete at sea level (e.g., Heraklion), we also need Evaporative loss there (in mm/day). The loops below determine a matrix of Evap values at 100 ALTITUDE increments up to 1000 m altitude, so we could use results for higher locations too:
For Crete:
For Haifa:
For altitudes of 0, 100, 200 ...1000 m above Heraklion and Haifa, we can now determine the E-P balance:
For Crete:
in mm
For Levant:
in mm
Net E-P balances (in mm) at sea level (k=0)
Control YEAR
Control SUMMER
Control WINTER
0 Event YEAR
0 Event SUMMER
0 Event WINTER
3 Event YEAR
3 Event SUMMER
3 Event WINTER
5 Event YEAR
5 Event SUMMER
5 Event WINTER
Because all events are in winter, the summer values don't change between different years.
Note that all values are in mm (not in mm/y), so annual values need to be obtained by addition of summer and winter values.
Now use Lapse Rates and surface temperatures to approximate snow-line altitudes in Crete and Haifa/Levant. This is important because key evidence for the RCC events consists of Little Ice Age reports of snow and frost. Also today, outbreak events are characterised by snow in unusual places (Athens, Istanbul).
Cretan snowline altitude
Levant snowline altitude
Lowest altitude of snowline during events
Xoplaki et al (2001 Clim. Cha.): "From the scattered data found for 1675-1715 and 1780-1830, the winter and spring climate in southern Balkans and the eastern Mediterranean, especially during the LMM, can be characterised as cooler and relatively rainier with a higher variability compared with the recent decades."
Snowline today in the western Taurus Mountains, Turkey is at about 3400-3700 m.
Cretan mountains reach up to almost 2500 m. Typically, Crete's three key mountain ranges reach above 2000 m. Event periods had persistent snow on the ranges, so snowline altitude 2000 m or lower.
At "cool" = 0, lowest snowline is at 2460 m. At "cool" = 1 it is at 2315 m. At "cool" = 2 it is at 2152 m. At "cool" = 3 it is at 1990 m. At "cool" = 4 it is at
Today, Thessaloniki air T range through year: 5 -26.5 deg C.
So to get
long-lasting frosts in winter in Thessaloniki
Xoplaki et al. (2001) write:
Note: Most extreme events :
If sustained/repeated in winter, this would give about 2.6 deg C SST change in winter relative to Control - again this is NOT impossible.
Use Rohling (1999 Paleoceanography) calculation for
Resistance coefficient Ratio = 1.0142 after Gonfiantini (1986), but some say this is too large (see Rohling 1999), so we use this to tune the 1-event "control" year rainfall deltas to about -6. I get close to this when I set Ratio = 1.0105 (for cool=6), or Ratio=1.0098 (for cool=10): we settle for 1.010
Aegean
Condensation fractionations for the O isotopes:
Initial rainfall
Isotope mass balancing between import, evap, and final
Calculate change (anomaly) in oxygen isotope values of Aegean surface water
Thus, the final
remaining
vapour
And ensure
Levantine
Initial rainfall
Now the entire
And ensure
Disable the functions when the code flags that "Sub" in the loops is empty
mean rainfall
Apparent interannual amount effect:
permil per 100 mm
Mean rainfall
Disable the functions when the code flags that "Sub" in the loops is empty
Apparent interannual amount effect:
permil per 100 mm
These runs use the cut-offs in thef GISP2 K+ ion series to identify periods with certain frequencies of northerly outbreaks. They allow for periods with fewer events than today by setting the control scenario to the 1-event years, and the 0-event years then are "more optimal" than today.
Categorising event frequencies from K+ scores:
50th percentile
80th percentile
95th percentile
Set counter value (h) for present-day value
All normalisations are OK as long as this value is equal to 1
Check that present-day is a 1-event year
First, a little
Set preferred altitude
Now make time-series using 0,1,3,5 to select the year characteristics as defined above.
Define annual values
Determine Haifa carbonate O-isotopes (using
22
Determine smoothed records to approximate recording in natural archives
Set smoothing factors (y):
Although these three anomalies are cumulative over the 3 years of each dt step, we are interested in their annual anomalies, which is what we'd record if we assume the system rests them each year.
Some records we want in anomalies, relative to present value:
carbonate O-isotope anomaly for 3k event re. to present:
Number of degree days record is multi-decadally smoothed to allow for longevity of snow fields (potential to exist through summers)
We can see that the annual freshwater budgets per single square meter are always negative in both locations. So to have sustained lakes, both locations need catchment areas that are larger than the areas of the lakes themselves:
Cretan
We grossly simplify by assuming that all P immediately soaks deep into the soil, and thus that in-catchment wE=0. wE only has value over the lake itself.
lake area:
catchment area
steady state:
during a "control" year (a 1-event year)
or override:
Lake-level change are in m/y
relationship used:
Lake-wall slope in degrees:
Initial lake depth: Set once.
in radians
Lake diameter (m):
diameter of lake
initial area of lake
area of lake+catchment
Levantine
We grossly simplify by assuming that all P immediately soaks deep into the soil, and thus that in-catchment wE=0. wE only has value over the lake itself.
lake area:
catchment area
during a "control" year (a 1-event year)
steady state:
or override:
Lake-level changes are in m/y
Lake-wall slope in degrees:
Initial lake depth: Set once.
in radians
Lake diameter (m):
diameter of lake
initial area of lake
area of lake+catchment
Cretan lake
Note: to determine dt
NB. All E and P values are given in mm, and we want answers for LL in m; hence divide by 1000
Repeat values:
set dt manually based on tdif determined above
Levantine lake
NB. All E and P values are given in mm, and we want answers for LL in m; hence divide by 1000
Repeat values:
Dead Sea specific case
For parameter settings, see Rohling (2013 QSR)
when LL>100 m (i.e., >-18m relative to present)
when LL<100m (i.e., =<-18m relative to present)
current evap reduction factor due to S; at near modern lake levels this is close to 0.5 (see Rohling 2013)
ratio of
0.025 assumes that lake will fill catchment area linearly with changes in height (1/40)
30
KEY RESULTS
Max. Aegean mixed layer cooling (deg. C):
Minimum snowline altitudes:
Apparent amount effects:
permil per 100 mm
me
Control year Precip values:
Mean annual O-isotope values for different types of year:
NB., "one" is the control scenario!
Set counters to pick specific columns
Vinc = 10 would be reasonable as Bora winds can be consistent at 17 m/s http://www.istrianet.org/istria/geosciences/meteorology/winds/bora-adriatic.htm
Then:
O isotope anomaly advected into region with incoming air during events (NB can be source-water effect due to NAtl 'hosing'):
Check for Thessaloniki frosts during the events: