SHARE
TWEET

Untitled

a guest Dec 13th, 2018 59 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy as np
  2. def K0_Weiss(S, TC):
  3.     TK=TC+273.15 #Temperatura de Celsius a Kelvin  
  4.     lnK0 = 9345.17/TK - 60.2409 + 23.3585 * np.log(TK/100) + S * (0.023517 - 0.00023656 * TK + 4.7036e-07 *
  5.         TK * TK)
  6.     K0=np.exp(lnK0)
  7.     return K0
  8.  
  9.  
  10. """
  11. Created on Mon Sep  3 11:51:00 2018
  12.  
  13. @author: fergomez
  14. Esta funcion calcula K1 y K2 del sistema CO2-H20 segun
  15. Roy et al., (1993), The dissociation constant of carbonic acid in sewater
  16. at salinities of 5 to 45 and temperatures of 0 to 45 desgrees Celsius.
  17. Marine Chemistry 44(2-4), 249-267.
  18. Unidades
  19. K1, K2=mol*kg-soln-1
  20. Ejemplo de uso:
  21. Input
  22. K1, K2=K1_K2_Roy(35, 25)
  23. print(K1, K2)
  24. Output:
  25. 1.3921075396202872e-06 1.1887254858040348e-09
  26. """
  27.  
  28.  
  29. def K1_K2_Roy(S, TC):
  30.     TK=TC+273.15 #Temperatura de Celsius a Kelvin    
  31.     #K1 usando Roy et al., 1993
  32.     tmp1 = 2.83655 - 2307.1266/TK - 1.5529413 * np.log(TK)
  33.     tmp2 = -(0.207608410 + 4.0484/TK) * np.sqrt(S)
  34.     tmp3 = 0.08468345 * S - 0.00654208 * S**1.5
  35.     tmp4 = np.log(1 - 0.001005 * S)
  36.     lnK1roy = tmp1 + tmp2 + tmp3 + tmp4
  37.     K1 = np.exp(lnK1roy)
  38.     #K2 usando Roy et al., 1993
  39.     tmp1 = -9.226508 - 3351.6106/TK - 0.2005743 * np.log(TK)
  40.     tmp2 = (-0.106901773 - 23.9722/TK) * np.sqrt(S)
  41.     tmp3 = 0.1130822 * S - 0.00846934 * S**1.5 + np.log(1 -
  42.         0.001005 * S)
  43.     lnK2roy = tmp1 + tmp2 + tmp3
  44.     K2 = np.exp(lnK2roy)
  45.     return K1, K2
  46.  
  47. def K1_K2_SBY(S):
  48.     pK1=6.1568-0.00352*S
  49.     pK2=8.5503-0.0080*S
  50.     K1=10**(-pK1)
  51.     K2=10**(-pK2)
  52.     return K1, K2
  53. """
  54. Created on Mon Sep  3 11:53:53 2018
  55.  
  56. @author: fergomez
  57. Esta funcion calcula K0 del sistema CO2-H20 segun recomienda Dickson et al 2007
  58. Guide for best practices for ocean CO2 measurements. PICES SP 3, 191 pp.
  59. Conforme a Mucci, A., 1983. The solubility of calcite and aragonite in seawater
  60. at various salinities and temperatures and one atmosphere of total pressure.
  61. American journal of Science, 283, 780-799.
  62.  
  63. Ejemplo de uso:
  64. Input
  65. Ksp_a=Kspa_Mucci(35, 25)
  66. print(Ksp_a)
  67. Output:
  68. 6.48175906801198e-07
  69. """
  70.  
  71. def Kspa_Mucci(S, TC):
  72.     TK=TC+273.15 #Temperatura de Celsius a Kelvin
  73.     tmp1 = -171.945 - 0.077993 * TK + 2903.293/TK + 71.595 * np.log10(TK)
  74.     tmp2 = +(-0.068393 + 0.0017276 * TK + 88.135/TK) * np.sqrt(S)
  75.     tmp3 = -0.10018 * S + 0.0059415 * S**1.5
  76.     log10Kspa = tmp1 + tmp2 + tmp3
  77.     Ksp_a=10**(log10Kspa)
  78.     return Ksp_a
  79.  
  80.  
  81. """
  82. Created on Mon Sep  3 11:54:42 2018
  83.  
  84. @author: fergomez
  85.  
  86. Esta funcion calcula K0 del sistema CO2-H20 segun recomienda Dickson et al 2007
  87. Guide for best practices for ocean CO2 measurements. PICES SP 3, 191 pp.
  88. Conforme a Mucci, A., 1983. The solubility of calcite and aragonite in seawater
  89. at various salinities and temperatures and one atmosphere of total pressure.
  90. American journal of Science, 283, 780-799.
  91.  
  92. Ejemplo de uso:
  93. Input
  94. Ksp_c=Kspc_Mucci(35, 25)
  95. print(Ksp_c)
  96. Output:
  97. 4.2723509278626e-07
  98. """
  99.  
  100. def Kspc_Mucci(S, TC):
  101.     TK=TC+273.15 #Temperatura de Celsius a Kelvin
  102.     tmp1 = -171.9065 - 0.077993 * TK + 2839.319/TK + 71.595 * np.log10(TK)
  103.     tmp2 = +(-0.77712 + 0.0028426 * TK + 178.34/TK) * np.sqrt(S)
  104.     tmp3 = -0.07711 * S + 0.0041249 * S**1.5
  105.     log10Kspc = tmp1 + tmp2 + tmp3
  106.     Ksp_c=10**(log10Kspc)
  107.     return Ksp_c
  108.  
  109. # -*- coding: utf-8 -*-
  110. """
  111. Created on Mon Sep  3 11:55:56 2018
  112. Calcula el factor de correccion de la presion para las constantes de
  113. equilibrio Ki usando Millero, F.J., (1995). Thermodynamics of the carbon
  114. dioxide system in the oceans. Geochimica et Cosmochimica Acta 59:661-677.
  115. @author: fergomez
  116.  
  117. Ejemplo de uso:
  118. A presion atmosferica Kspa=4.27 e-07
  119. Input (para aragonita):
  120. fcorr=P_corr_Millero(-45.96, 0.5304, 0.0, -11.76e-03, 0.3692e-03, 0.0, 25, 300)
  121. print(fcorr)
  122. Output:
  123. 1.4787577790538065
  124.  
  125. entonces Kspa(presion atmosferica)*fcorr=Kspa_corregido
  126. 4.27 e-07*1.4787577790538065= 9.58 e-07
  127. """
  128.  
  129.  
  130. #pagina 186 Andreas Hoffman PhD Thesis
  131. def P_corr_Millero(a0, a1, a2, b0, b1, b2, TC, P):
  132.     TK=TC+273.15 #Temperatura de Celsius a Kelvin
  133.     R=83.131
  134.     delta_V= a0+a1*TC+a2*TC*TC
  135.     delta_k= (b0+b1*TC+b2*TC*TC)
  136.     lnKp=-(delta_V/(R * TK))*P+0.5*(delta_k/(R*TK))*P*P
  137.     fcorr=np.exp(lnKp)
  138.     return fcorr
  139.  
  140. #Llamamos las librerias que utilizaremos
  141. import numpy as np
  142. #import matplotlib.pyplot as plt
  143. import math
  144.  
  145. from numpy import exp,arange
  146. from pylab import meshgrid,cm,imshow,contour, clabel,colorbar,axis,title,show, contourf
  147.  
  148. #from constantes import K0_Weiss, K1_K2_Roy, Kspc_Mucci, Kspa_Mucci, P_corr_Millero
  149. import matplotlib.pyplot as plt
  150.  
  151. #============================================================================
  152. #Inicio de seccion de INPUT==================================================
  153. #============================================================================
  154.  
  155. S=6.88 #ppt
  156. TC=6.42 #Celsius
  157. P=1.0 #atmosferas
  158.  
  159.  
  160.  
  161. K0=K0_Weiss(S, TC)
  162. K1, K2=K1_K2_Roy(S, TC)
  163. Ksp_c=Kspc_Mucci(S, TC)
  164. Ksp_a=Kspa_Mucci(S, TC)
  165.  
  166. #============================================================================
  167. #Fin de seccion de INPUT======================================================
  168. #============================================================================
  169.  
  170. #=============================================================================
  171. #Estimacion de Calcio de Tyrrell et al., 2008 y valores de ejemplo agua de mar
  172. #============================================================================
  173. #Ca=0.0103 #valores agua de mar
  174. Ca=(0.331*S+0.392)*1e-03 #Mar Baltico
  175. #Ca=(0.375*S+0.0368)*1e-03 # Bothnian Bay (parte norte del mar Baltico)
  176. #=============================================================================
  177.  
  178. #def co2_t(pco2, K0):
  179. #    co2=pco2*K0
  180. #    return co2
  181. #=============================================================================
  182. #vco2_t = np.vectorize(co2_t)
  183. #Pasamos de pco2 a co2
  184. #co2=co2_t(pco2, K0)
  185.  
  186.  
  187. #=============================================================================
  188. #=============================================================================
  189. #Calculo del factor de correccion por presion
  190. #=============================================================================
  191. fcorr_K1=P_corr_Millero(-25.50, 0.1271, 0.0, -3.08e-03, 0.0877e-03, 0.0, TC, P)
  192. fcorr_K2=P_corr_Millero(-15.82, -0.219, 0.0, 1.136e-03, -0.1475e-03, 0.0, TC, P)
  193. fcorr_Kspc=P_corr_Millero(-48.76, 0.5304, 0.0, -11.76e-03, 0.3692e-03, 0.0, TC, P)
  194. fcorr_Kspa=P_corr_Millero(-45.96, 0.5304, 0.0, -11.76e-03, 0.3692e-03, 0.0, TC, P)
  195. #=============================================================================
  196.  
  197. #=============================================================================
  198. #Correccion de ctes de equilibrio por presion
  199. #=============================================================================
  200. K1=K1*fcorr_K1
  201. K2=K2*fcorr_K2
  202. Ksp_c=Ksp_c*fcorr_Kspc
  203. Ksp_a=Ksp_a*fcorr_Kspa
  204.  
  205. #=============================================================================
  206. #Funcion que realiza calculos del sistema carbonato y omega
  207. #=============================================================================
  208.  
  209. def carbEq(co2, dic, Ca):
  210.     #-----------------------------------------------
  211.     # Resolvemos para obtener H+ (cf. a Zeebe and Wolf-Gladrow, 2000)
  212.     a1=co2-dic
  213.     a2=K1*co2
  214.     a3=K1*K2*co2
  215.  
  216.     p= [a1, a2, a3]
  217.     r = np.roots(p)
  218.     h= max(np.real(r)) # Esto es para seleccionar la raiz real mas grande
  219.     #
  220.     # Calculamos HCO3, CO3 and CO2aq, usando DIC, AlK y H+
  221.     hco3=dic/(1+h/K1+K2/h)
  222.     co3=dic/(1+h/K2+h*h/(K1*K2))
  223.     #co2=dic/(1+K1/h+K1*K2/(h*h))
  224.     fco2=co2 / K0
  225.     pH=-math.log10(h)
  226.     alk=2*co3+hco3
  227.  
  228.     #Saturacion de Calcita y Aragonita
  229.     omega_ar=Ca*co3/Ksp_a
  230.     omega_cal=Ca*co3/Ksp_c
  231.  
  232.     return fco2, pH, co2, hco3, co3, h, dic, alk, omega_ar, omega_cal
  233.  
  234. vcarbEq = np.vectorize(carbEq)
  235. #=============================================================================
  236. #Uso de la funcion - ejemplo
  237. #=============================================================================
  238. #co2=co2*1000000
  239. #dic=dic*1000000
  240.  
  241. pco2=np.linspace(0.000090, 0.001100, 100)
  242. co2=pco2*K0
  243. dic=np.linspace(0.001526, 0.002100, 100) #moles (para pasar de micromoles a moles, 1587*e-06)
  244.  
  245. #pco2=pco2*1000
  246. #dic=dic*1000
  247. co2, dic=meshgrid(co2*1000000, dic*1000000)
  248.  
  249. fco2, pH, co2, hco3, co3, h, dic, alk, omega_ar, omega_cal=vcarbEq(co2=co2, dic=dic, Ca=Ca)
  250.  
  251.  
  252. #Ejemplo con valores de referencia de uso para agua de mar normal
  253. #fco2, pH, co2, hco3, co3, h, dic, alk, omega_ar, omega_cal=carbEq(co2=0.00001032997, dic=0.002108, Ca=Ca)
  254. #para este usar calcio=0.0103
  255. #=============================================================================
  256. #Grafico de los resultados
  257. #=============================================================================
  258. print(co2)
  259. print('====================================')
  260. print(dic)
  261. print('====================================')
  262. print(pH)
  263. #print(omega_cal)
  264. #fig, ax = plt.subplots()
  265. im = imshow(pH, cmap=cm.RdBu) # dibujo la funcion
  266. plt.xlabel('CO2')
  267. plt.ylabel('DIC')
  268. plt.ylim(0, 100)
  269.  
  270. plt.scatter(40, 40, color='k')
  271. #ax.set_xlim(0.00000525, 0.000065)
  272. #ax.set_ylim(0.001, 0.0021)
  273. #plt.xlim(min(co2), max(co2))
  274. # agrego lineas de contorno y rotulos
  275. cset = contour(pH, arange(6.0,9.0,0.2),linewidths=2, cmap=cm.Set2)
  276. clabel(cset,inline=True,fmt='%1.1f',fontsize=10)
  277. colorbar(im) # agrego la barra de colores al costado
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top