PTPath + Isentropic

[robert]

Hello,

I notice when I run adiabat in isentropic mode along a constant dP the temperature drops gradually from it's starting point, which I assume is due to thermal energy being lost at the onset of melting and phase changes. I'm wondering if there's some way to combine this temperature decline with a predetermined PTPath.

Essentially, I have a PTPath of the upper mantle with temperatures defined by purely conductive cooling. If I run adiabat along this path from depth to surface I'm afraid I'm neglecting the effective temperature drop caused by phase transitions (etc). What I would like is a way to incorporate this temperature drop with my conductive cooling temperature model somehow, but in order to do so it seems like what I would have to do is alternate between PTPath and isentropic modes, every other step.

If I may further explain,

say my PT values for PTPath are:

26000 1350
25000 1348
24000 1345
...

and the results in isentropic mode, given the same initial T:

26000 1350
25000 1349.5
24000 1348.8
...
.......which decrease only slightly in temperature, and only because melting is occurring (and at significant phase boundaries). I think?

If anyone has a scheme or could possibly make this possible in a future version, or if it's currently possible now and I'm missing something, I'd like to see the thermal effects of melting incorporated WITH the PTPath, making it a conductive AND thermodynamic temperature profile. In other words, allowing adiabat to fudge my PTPath values in order to incorporate the additional thermal effects.

My guess is it'd be something like:

26000 1350
25000 1347.5
24000 1343.8

..but not exactly a direct superposition of the different outputs since I assume melting is somewhat temperature dependent  ;)

Currently my idea would be to run in isentropic, but adjust the starting temperature before each iteration based on what my conductive cooling model predicted. Ideas?

Thanks in advance!

[Paula]

Hi Robert,

Just a quick reply for now:

Some of the variation in temperature in isentropic mode will be due to melting or other phase changes but much of it (~ 2 ^oC per kbar for typical peridotite) will be simply adiabatic decompression. Your conductive cooling P-T paths are presumably for an incompressible solid(?). You'll be able to distinguish between temperature changes due to decompression or melting to first order by supressing liquid, though eliminating sub-solidus reactions is a little more difficult.

It should be possible to account for the latent heat of melting (or other reactions). The public version of pHMELTS doesn't have a direct way of doing that, at least not that I can think of offhand. Adjusting the temperature at each iteration is unlikely to work as either (a) if it's the first isentropic calculation and it will use the current temperature to set the reference entropy or (b) if it's a later calculation and the reference entropy is already set then it will use it so your temperature adjustment will have no effect.

More later,
Paula

[robert]

Right- essentially it's an incompressible halfspace model, no melt, constant viscosity, not boundaries, as simple as possible. I feel the melt-related temperature drop over the range of depths I'm attempting to model is significant (drops 100+ C from a range of 26 kbar or so) and would like to somehow incorporate that into my PTPath.

To my knowledge there isn't a setting for this on the vanilla public version, and it does seem that any roundabout way I can do this will be flawed, so if there's some way to allow for PTPath + appropriate temperature fudging I would be very grateful to see it in the next version. Thanks!

[Paula]

If you tell me how you want the 'appropriate temperature fudging' to be calculated (or input?) then I can look into implementing it. But I think you need to be a bit clearer about what you are trying to model...  Presumably at the axis (time=0) you impose a uniform temperature with depth (plus a top boundary condition) and then you calculate the conductive profile at a given time, or rather distance from the axis, either analytically or from a parameterisation. Does the model have a smooth thermal boundary layer?

I'm guessing that what you'd like to get out is a temperature profile for each time step as if the initial temperature profile had accounted for melting? If you are not extracting the melt then eventually you will get back to something close to your basic model as the latent heat of crystallisation will balance the budget, so to speak, in the far-field. Adiabat_1ph can't tell you how the temperature 'anomaly' due to having melting at the axis will evolve in the short term; wouldn't that require solving the (2-D?) heat equation numerically? If you're extracting melts and adding crust then it's even more complicated as the heat sinks and sources are no longer in the same place.

You can probably estimate what would happen by allowing the temperature difference between two Adiabat_1ph calculated on-axis profiles with and without liquid suppressed to dissipate over a characteristic time scale, whilst always fixing the base of the plate (sensu lato) to the 'with-melt' calculated temperature. But such an approximation is easier done outside of Adiabat_1ph with some sort of script. It might make things easier if your thermal model had an adiabatic temperature gradient below the plate (e.g. see Ritzwoller et al.)? On the other hand, having a thermal boundary layer (or 'transition layer' as they call it) could make the calculation a little harder as the point at which melting switches off may be less well defined.

If I'm totally missing what you are suggesting then please let me know in a bit more detail.

Paula

[robert]

I know I'm explaining this poorly, I'm obviously no petrologist. But let me try again (and again if needed).

I have a vertical, 1D column P,T profile of the upper mantle, derived from halfspace model. Any change in temperature is due to conductive cooling alone, so there aren't any definite boundary layers to speak of. I want to use this data as a ptpath from depth to surface in order to determine melt quantities and compositions at each 1km depth interval, assuming fractional melting.

Now, when melting occurs at each depth, temperature at that depth should decrease a little (again, assuming fractional melting). And I can see by running in isentropic mode that indeed it does, and it's not negligible for the wide range I'm going through. 100C per 25kbar, thereabouts.

My end goal is not to come away with a valid T profile as much as an accurate representation of what's being melted. So, if I force each step to be strict to my "simple", conductive-only ptpath, the melt data will be off because the temperature [at all depths other than my starting point] will be higher than it should, because all the pervious melting and phase transitions will have taken up energy.

So, in a perfect world (for me anyway), adiabat can "officially" run in ptpath mode but be aware of how much temperature is being lost at each step and apply that change in temperature to the next step, as it does in isentropic mode.

I think you know what I'm getting at since you wrote I'm guessing that what you'd like to get out is a temperature profile for each time step as if the initial temperature profile had accounted for melting? because If I had *that* ptpath I'd be set. Sort of like running in isentropic mode, calculating the dT for each pressure interval, then retroactively applying the resulting isentropic dT profile to my ptpath and running from that.

I'm not really thinking about this in terms of time steps really, nor am I trying to sum two dimensionally and/or for a crustal thickness or anything. I just need f and melt compositions at each point on my ptpath that I'm fractionally extracting for use in another code.

Thanks for your patience!
Robert 

[Paula]

I think you know what I'm getting at since you wrote I'm guessing that what you'd like to get out is a temperature profile for each time step as if the initial temperature profile had accounted for melting? because If I had *that* ptpath I'd be set. Sort of like running in isentropic mode, calculating the dT for each pressure interval, then retroactively applying the resulting isentropic dT profile to my ptpath and running from that.

No, I don't know what you're getting at! I though I did previously but, given that you now say you are not really interested in the temperature profile, I can't work out if your comment is a typo or sarcasm. If you want melt fractions and compositions then you would probably not be 'set' with that PT path as it will be too cold for anything to melt (unless the lithosphere had been enriched metasomatically but that's a different story). As you mentioned 25-26 kbar a couple of times, I am guessing that your 1-D model is for a plate thickness of ~85 km, corresponding to an age of ~70 Ma. It is a snapshot in time or space, depending on your frame of reference, of the half-space cooling model but any melting occurs elsewhere or earlier and in a different thermal regime:

This figure is adapted from Figure 1 in Asimow & Langmuir for two different potential temperatures. I know things are not quite to scale but it illustrates the point. If the blue lines represent your 1-D model then, if I understand correctly now, what you effectively need to know is the residual composition for that vertical column. So what you want, I think, is the aggregate melt fractions and compositions along the top surface of the melting region i.e. as the streamlines leave the melting region.

The figure shows a triangular melting region but for a simple model in which the streamlines turn a sharp corner it makes no difference in the far-field whether the melting region is perfectly triangular or concave up etc.. Melting will not occur spontaneously in regions where material is moving horizontally, even without conductive cooling. At the time when the melt is actually extracted the mantle hasn't had time to cool anyway so you should be able to calculate what you need with a 1-D isentropic fractional (or continuous) melting calculation i.e. along the red line. You can use the standalone or internal version of cubicquad2 to find the 1-D aggregate melt fraction and composition at each pressure. (If there is a thermal boundary layer below the plate, so that there is not an abrupt change from adiabatic decompression to horizontal motion, then conductive cooling will have an effect on melting. But that is more complicated than what you are trying to capture in your model and would be beyond the scope of Adiabat_1ph, as discussed before.)

If you needed to know how the assemblage and phase compositions of the residual mantle column change as it re-equilibrates during lithospheric cooling then that is, indeed, not something included in the public version of Adiabat_1ph at the moment. It is possible, but would be time consuming to get into a reasonably user-friendly form for public use (and I know you have a fairly low opinion of the current user-friendliness, which I'm still working on...). If you really only want melt fractions and compositions then this would be superfluous to your needs anyway.

Hope that answers your question but, if not, please let me know.
Paula

[robert]

Hi Paula, me again.

No sarcasm, but I wouldn't be surprised if I blatantly contradicted myself.

To reiterate (for myself, I've been away for some time), I'm interested in the melt produced at 1 km intervals along a vertical path below a MOR, at all different distances from the axis (/times), some of which won't necessarily be erupted. I have many PTPath files for doing this, and it works pretty well. BUT- these files are derived from a strictly conductive cooling model and thus are inherently inaccurate in regions where melt exists, as the melting process will lower temperatures. If I run in isentropic mode starting at the initial P,T from my ptpath files, I notice that the code will lower the T values accordingly at the onset and continuation of melt. What I'm looking for is a way to run in PTPath mode, at my rigid predefined conductive cooling model, but incorporating the dip in T caused by melting. So the T values in my ptpath files are not exactly set in stone so to speak, they can be pre-altered before the next calculation in the same way as in isentropic mode. I think the numbers in my initial post may still help explain it.

With the figure you posted, I would agree that just running in isentropic mode should be sufficient, as temperatures along the red arrows should remain constant and cool only as melt is produced. But in the model we're using this isn't necessarily the case, as conductive cooling is a huge factor. Here's an edited version of that file if it helps:



In the center square I've drawn a typical vertical temperature profile in our model. At point "A" it's essentially isentropic mode, at point "B" is a transition to a region of both conductive cooling (our ptpath) AND "melting cooling" (isentropic mode), and by point "C" there should be no more melting occuring.

So, I guess we're shooting for more accuracy than can be afforded by just an isentropic model, and more still than the PTPath mode because then we're neglecting the temperature decrease melting causes. At this point all I can think to do is run one step in isentropic, then subtract the dT from the temperature in the next ptpath file step, run in PTPath mode one step, then repeat until we're at the surface.

Also, I tested out the nocheck/celsius versions of adiabat_1ph and they seem to work great! Thanks!

[Paula]

Hi Robert,

Thanks very much for the testing those two things. I'm just tying up the loose ends to get v1.9 up - should be done by some time tomorrow. I like the picture editting! I'm still a bit unsure about some stuff though. Please could you do me a favour and either post or e-mail me two or more example PT files from you model: one from near, but not quite on, the ridge axis and one from further away. Make sure it's the whole vertical section including the 'isentropic' (really isothermal...?) bit.

Thanks,
Paula

[Paula]

Hi Robert,

Thanks for the PT profiles and information. I have tried using them in Adiabat_1ph with a pared down version of the Workmann & Hart bulk composition (removed NiO and other oxides that are not calibrated in pMELTS; also removed K2O as there's no solids for that to go into at high pressure). The thick black line below is the result of an isentropic calculation in Adiabat_1ph and the yellow star is the point at which melting begins (when there is no K2O in the bulk composition). So the figure illustrates what I was saying initially about most of the cooling being due to adiabatic decompression, rather than melting. It is not possible to 'turn off' this effect. As described in Mark Ghiorso's MELTS documentation, only certain constraints are possible in the thermodynamic calculations and S and V are not amongst the allowed combinations. (Adiabat_1ph does have a ADIABAT_IMPOSE_FO2 enviroment variable that allows you to approximate the '+ fixed fO₂' option for all boxes in the table that have a 'yes', but that's an aside.)

So, I would suggest tailoring the PT paths that your cooling model calculates instead. The PT paths you provided are the dotted lines and the solid coloured lines are schematic profiles for the cooling model corrected for the adiabatic gradient. I did this is a graphics program by shearing and scaling but it's not hard to see how you could code up something that would give a pretty good approximation of the effect of (pMELTS-calculated) adiabatic decompression on the cooling profile. You could, if you want, include the change in temperature due to melting to make a composite isentropic and conductive cooling path like the dashed black line. As I've said before, this is really something that should be done outside of Adiabat_1ph (using the 1-D isentropic path as input). If I was doing the bit of research you are trying that is how I would do it, either within the source code for the thermal model or by writing a script to process the thermal model output before passing it to Adiabat_1ph. As discussed earlier in the thread there are some assumptions (e.g. ignoring any latent heat of crystallisation later on, rate of diffusion of the thermal 'anomaly') but it would be a start.

For the bulk composition I used, no melting occurs in the off-axis profile. When there is K2O is included in the bulk composition, there is a tiny amount of melt present almost everywhere (both profiles) as there are few solids in the MELTS database, other than feldspar, that incorporate K. But the amount of melt does not vary much so the latent heat of melting would be negligible and, besides, in the cold profile the amount of melt is well below the default threshold for extraction. Incidently, 1300 ^oC is rather low for 'normal' mantle potential temperature in the pMELTS world (due to systematic offsets e.g. see Hirschmann et al 1998,  Ghiorso et al 2002, Asimow et al 2001; full references are in the documentation). You can get a more appropriate estimate by doing a ordinary 1-D isentropic calculation with the cubicquad2 integration (option 12) to get the thickness of the 2-D aggregated crust. With a bit of trial and error you can choose the temperature that gives typical crustal thicknesses for the area you are studying. Depending on the composition (which oxides are excluded, whether you include water as a trace element pHMELTS-style etc.) a value closer to 1400 ^oC may be better.

If you wanted to model curved streamlines like the ones you drew then you need to rearrange your PTpath files (not that I am recommending that you do, by the way). The ptpath_file is supposed to be just that: a path. For batch melting it doesn't matter but fractional melting is irreversible and so the path should (normally) trace a packet of mantle along the streamline it follows. If the mantle follows the blue streamlines, then the bulk compoistion of the piece of mantle at 1 is not related to that at 2 by any simple melting relationship. How is Adiabat_1ph supposed to know the history of the mantle at each depth in the 1-D profile?

On the other hand, you could assume uniform upwelling and have the streamlines turn a sharp corner at the transition from adiabatic upwelling to conductive cooling. Then at a given depth in the melting region all the residual material has the same bulk composition (RHS figure) and a vertical profile can be used. The 1-D isentropic calculation would make up the upwelling part of the profile, similar to what was proposed above. In reality, neither 1 nor 2 is the same as the original bulk composition, as both will have been depleted by melting on-axis. So your current off-axis profiles may overestimate the amount of melting not just because the temperature has not been adjusted for adiabatic decompression and / or melting but also because the composition will be too fertile. Now, if you are only interested in the amount and composition of melt generated then you don't even need to do the calculation in the conductively cooled part, so the bulk composition is not an issue. You don't need to assume a triangular melting region like the figures here - you could use your cooling model to calculate the plate thickness and stop the calculation at the base of the plate each time.

If, hypothetically, you were interested in the mineral assemblage or compositions in the residue after cooling then that would require adding something like a 'PTPath + Isentropic' mode to Adiabat_1ph. As I alluded to before, I have done a calculation like that once with a Caltech colleague but the current procedure is user-unfriendly and labour intensive (and we only did one profile). Although I can now see a possibly better way to implement it, you have to understand that it would only be able to adjust the residual composition. The temperature of the PT path would still need to be adjusted for adiabatic decompression and melting, if you want. But you would have to do that yourself as set out above. As you have described it, the 'PTPath + Isentropic' mode would not help you with your current problem but if you, or anyone else, were interested in calculating the off-axis residual mantle column more easily then I would consider adding this functionality in future .

Paula