Simulação de perfis de espectrometria de raios gama naturais em meios aleatórios

Simulation of natural Gamma-ray spectrometry logs in random media

Simulation of natural Gamma-ray spectrometry logs in random media

J. da Conceição da Silva

Universidade Federal do Rio de Janeiro

Departamento de Geologia - Setor de Geologia de Engenharia e Ambiental

Instituto de Geociências - Prédio do CCMN, BL - J0

Av. Brigadeiro Trompowiski, s/n, Ilha do Fundão

CEP.: 21.949 - 900, Rio de Janeiro, RJ.

E-mail: [e-mail]

In this work a theoretical model has been developed to simulate the natural gamma-ray spectrometry responses in random media. The transport of gamma-rays is simulated by the finite-element method (FEM) for the scalar flux derived from the P₁ approximation of the Boltzmann transport equation for radially symmetric media. The spectral response concerned with the random distributions of potassium-40, uranium-238 and thorium-232 enables us to determine the concentrations by weight of these nuclides using the Monte Carlo method. The physical parameters of the formations, i.e., diffusion coefficient and total and scattering cross-sections, are determinated by Klein-Nishina distributions and multiple linear regression of the experimental data. This method is used to generate synthetic well logs and also shows the usefulness of the FEM technique in extracting important information from the nuclear logs run in a complex geological environment.

Key words: Well-logging; Gamma-ray spectrometry; Random sources.

Simulação de perfis de espectrometria de raios gama naturais em meios aleatórios - Neste trabalho, é desenvolvido um modelo teórico para simular o perfil de espectrometria de raios gama naturais em meios aleatórios. O transporte da radiação gama é simulado através do método de elementos finitos (FEM) que calcula o fluxo escalar proveniente da aproximação P₁ da equação de transporte de Boltzmann para meios radialmente simétricos. A resposta espectral baseada na distribuição randômica de potássio-40, urânio-238 e tório-232 permite determinar as concentrações em peso destes nuclídeos conforme sugerido pelo método Monte Carlo. Os parâmetros físicos das formações, isto é, o coeficiente de difusão e as seções de choque total e de espalhamento, são determinados através da distribuição de Klein-Nishina e de regressão linear múltipla de dados experimentais. O método é então aplicado na geração de perfis sintéticos de poço e mostra a utilidade da técnica de elementos finitos na predição de informações importantes contidas nos perfis nucleares corridos em ambientes litológicos complexos.

Palavras-chave: Perfilagem de poços; Espectrometria de raios gama; Fontes aleatórias.

INTRODUÇÃO

A identificação de formações geológicas em subsuperfície constitui-se hoje em um dos mais intrigantes desafios impostos pelas modernas técnicas de interpretação geral de dados geofísicos e avaliação de formações. Uma forma corrente de obter parte destes dados geofísicos consiste em determinar os conteúdos de potássio, urânio e tório, ou, em alguns casos, o conteúdo total destes radioelementos, fazendo baixar um detector de radiação gama em um poço que intercepta os horizontes a serem investigados e medir a sua radioatividade natural. Além de ser uma técnica simples de perfilagem, este método permite contornar sérios problemas de ambigüidade apresentados por outros métodos geofísicos (Meyers, 1992). A eficiência dessa técnica está no fato de que as distribuições espectrais de K-40, U-238 e Th-232 estão associadas a uma litologia específica (Wahl, 1983).

Neste trabalho são simulados perfis sintéticos de espectrometria de raios gama naturais utilizando como fontes os radioisótopos citados acima, porém distribuídos aleatoriamente ao longo das camadas que constituem modelos geológicos complexos. Para tal, o transporte da radiação é calculado através do algoritmo de elementos finitos (FEM), na forma de fluxo escalar 2 1/2-D, proveniente da solução numérica da aproximação de difusão para multigrupos da equação de transporte de Boltzmann, no espaço de fase (Nowosad, 1978). Esta é conhecida como aproximação P₁, onde a variável direção e a fonte externa são expandidas em harmônicos esféricos através dos polinômios ortogonais de Legendre (Etherington, 1958).

As fontes isotópicas são simuladas por um processo simples de geração de fontes aleatórias baseadas no método Monte Carlo (Rubinstein, 1981). Neste método, os parâmetros físicos da formação são estabelecidos através de conceitos fundamentais de física nuclear aplicada e quântica não-relativística, de acordo com a distribuição de Klein-Nishina (Davisson & Evans, 1952). Mas, as concentrações em massa de K-40, U-238 e Th-232, em vez de apresentarem valores constantes para cada litologia, assumem valores distribuídos aleatoriamente dentro dos intervalos de variação observados experimentalmente para cada mineral presente nas rochas (Wahl, 1983). A seguir, a metodologia proposta é aplicada na elaboração de perfis sintéticos de poços com a finalidade de demonstrar a multiplicidade de aplicações deste importante sistema de perfilagem na avaliação de formações.

TRANSPORTE DE RADIAÇÃO

A equação de transporte de Boltzmann monoenergética, no estado estacionário e espaço de fase, para um sistema heterogêneo sem multiplicação, é dada por

(1)

onde W é o vetor unitário direção, r é o vetor espacial, S(r, W) é a fonte externa e Y (r, W ) é o fluxo angular (Barry & Pollard, 1982). A seção de choque total é dada por , onde N_a(r) e N_s(r) são as densidades de absorvedores e espalhadores presentes na formação, e s_a (W) e s_c (W) são as seções de choque microscópicas de absorção total e espalhamento Compton, respectivamente. A seção de choque diferencial microscópica, s_s (W.W’), define a probabilidade de espalhamento de fótons para a direção W a partir da direção W’.

A fim de adequar-se ao algoritmo FEM (Finite Element Method), a Eq. (1) deve ser reduzida pela aproximação P₁. Isto é feito expandindo as distribuições da fonte e do fluxo angular em harmônicos esféricos do cosseno do ângulo q (entre a direção W e a direção principal r) e do cosseno do ângulo a (entre as direções W’ e W (Etherington, 1958)). Após multiplicar P₁ por cada harmônico esférico P_n, integrando de -1 a +1, obtém-se uma série de n + 1 equações diferenciais de 1ª- ordem, de acordo com as propriedades de ortogonalidade dos polinômios de Legendre (Nowosad, 1978).

Na aproximação P₁ é assumido que, nas expansões, a contribuição dos termos após os dois primeiros são desprezíveis (Etherington, 1958). Estas considerações significam que o meio é fracamente absorvente e que os harmônicos de ordem superior expressam efeitos transientes próximos das fronteiras. Assim, a aproximação P₁ exige que a região em estudo esteja vários livres percursos médios distantes de quaisquer fronteiras.

Em geral, a aproximação P₁ toma a seguinte forma:

(2)

onde S₀ (r) é a densidade total de fótons da mesma energia que os de F₀, produzidos por fontes externas. O coeficiente de difusão é

(3)

onde s₁ é o segundo termo (i = 1) da expansão de s_s (W.W’) através de polinômios de Legendre P₁ (cos a).

Para o caso geral de difusão para multigrupos, supõe-se que as energias dos fótons estão distribuídas em um número finito de grupos, dentro dos quais as propriedades físicas do meio são constantes. Para descrever esta aproximação, o intervalo de energia [E₀,E_G] é dividido em G subintervalos I_g=[E_g-1,E_g],g = 1,2, ,G tal que, E₀ > E₁ > E_G com |E_g-1,E_g| ® 0 quando G ® ¥. A equação para multigrupos determina aproximações F_g (r), para o fluxo escalar exato integrado em cada intervalo de energia I_g, e assim a Eq. (2) torna-se

(4)

onde o último termo da Eq. (4) explica os termos da integral envolvendo Fg (r) que tinham sido incorporados á fonte aleatória S_g (r). Introduzindo agora uma seção de choque total preliminar vê-se claramente que

(5)

onde b_gg’ (r) é uma medida do número mais provável de fótons de energia E_g’ passando para a energia E_g, por unidade de massa, tempo e volume no ponto r. Para que o sistema de Eq. (4) não se decomponha em dois grupos, um dos quais não atua no outro, devemos supor b_gg’ (r) irredutível, responsável pelo acoplamento entre os "G" grupos (Nowosad, 1978).

Agora, introduzindo o operador linear como sendo os elementos da matriz diagonal L, considerando que F,s e B são os vetores fluxo escalar, fonte aleatória externa e matriz de transição, respectivamente, a equação de difusão é reduzida a um sistema de equações diferenciais na forma LF= BF+ s.

A equação final, para meios radialmente simétricos (Fig. 1),

Figura 1
- Modelo geológico no sistema de coordenadas cilíndricas.

Figure 1 - Geological model in cylindrical coordinate system.

assume a seguinte forma:

(6)

onde

,

, e

. Esta equação é resolvida por elementos finitos satisfazendo as seguintes condições de contorno: a simetria axial implica na condição de fronteira

(7)

a grandes distâncias da fonte, o fluxo deve aproximar-se de zero, isto é

(8)

e, finalmente, nas interfaces internas devem ser observadas as condições de continuidade do fluxo e da corrente. Uma completa descrição do método de elementos finitos pode ser encontrado em vasta literatura ( Oden & Reddy, 1976; Davies, 1980, entre outros).

PARÂMETROS FÍSICOS DA FORMAÇÃO

O coeficiente de difusão é determinado inserindo na Eq. (3) o valor de s1. Adotando as relações P₁ (m₀) = m₀ e m₀ = cosa , este parâmetro será

(9)

cuja solução requer apenas o conhecimento da seção de choque microscópica de espalhamento, s_s (cosa), o que pode ser obtido através da distribuição de Klein-Nishina (Davisson & Evans, 1952). Já o cálculo de b_gg’ (r,z) é feito integrando sob a curva de s_s^(g’) (q), entre os limites angulares [q_g’⁽¹⁾, q_g’⁽²⁾].

O esquema da Fig. 2 será tomado como base para o entendimento das diversas possibilidades de espalhamento.

Figura 2
- Esquema para o cálculo da probabilidade de transição de fótons entre os grupos de energia.

Figure 2 - Schematic diagram to calculate the probability of photon transition between energy groups.

Inicialmente, é analisado o caso em que o espalhamento se dá a um pequeno ângulo após o qual o fóton permanece no seu grupo de origem, isto é, g’ = g. Neste caso, q_g’⁽¹⁾ = 0. Já no caso mais geral onde g’ ¹ g, a fronteira superior do grupo g pode ser identificada calculando

(10)

onde mc² = 510,8 keV é a energia de repouso do elétron.

Se |cos q_g’⁽¹⁾| > 1, a probabilidade de espalhamento entre a energia do grupo anterior E_g’ e a energia E_g⁽¹⁾ relativa à fronteira superior do grupo g é zero, e assim, de b_gg’ (r, z) = 0. Quando |cos q_g’⁽¹⁾| < 1, deve-se pesquisar se o espalhamento se extende até o limite inferior do grupo g ou então identificar o ponto E_g’⁽²⁾ onde |cos q_g’⁽²⁾| £ 1. Pode-se determinar, analogamente,

(11)

Uma incoerência ocorre quando E_g⁽²⁾ > E_g’. Neste caso faz-se q_g’⁽²⁾ = p, valor este que se repetirá sempre que g’ for o último grupo, ou seja, aquele de menor energia, independentemente de ser |cos q_g’⁽²⁾| > 1 ou não. Enfim, de posse destes limites angulares, a seção de choque diferencial será dada por

(12)

onde s_s^(g’) (q) é a seção correspondente, obviamente, à energia do grupo g’.

A seção de choque total é composta por uma parcela que corresponde ao espalhamento Compton e outra devido a absorção de parte da energia do fóton nos processos de recúo do elétron Compton, absorção fotoelétrica e formação de pares.

O parâmetro é, então, a seção de choque macroscópica parcial para absorção de fótons, onde _as_s, t e k são as seções de choque microscópicas devidos aos processos de absorção citados acima, e N (r, z) é a densidade eletrônica da formação. A seção de choque a ss é calculada através da distribuição de Klein-Nishina, enquanto t e k são determinados por regressão linear múltipla dos dados de Davisson & Evans (1952), cujas equações são:

e (13)

Nestas expressões, válidas somente para 1£ Z_eff£ 20 , Eg é a energia dos fótons em keV, Z_eff é o número atômico efetivo da formação, definida de acordo com Willie (1984), e a e b são parâmetros ajustáveis em função da energia do fóton e do número atômico efetivo da formação. A Fig. 3 mostra uma boa concordância entre os dados experimentais e os resultados de t e k interpolados com a Eq.(13), para Z_eff=12, a=12,1 e b=1,06´10^-17 . A Fig. 3 mostra ainda a seção de choque microscópica total Compton de absorção e espalhamento intrínseco, s_C=_as_S+s_S, da distribuição de Klein-Nishina comparadas com os dados de Davisson & Evans (1952).

Figura 3
- Seções de choque microscópicas de absorção fotoelétrica (t), espalhamento Compton (s) e produção de pares (k).

Figure 3 - Photoelectric absorption (t), Compton scattering (s) and pair production (k) microscopic cross-sections.

A Tab. 1 mostra alguns dos valores numéricos para as seções de choque. Verifica-se que as seções de choque oriundas da regressão linear múltipla (RL) e da equação de Klein-Nishina (KL) não diferem mais que 3,9% dos valores (DE) determinados experimentalmente por Davisson & Evans (1952), o que confere a elas confiabilidade aceitável. Evidentemente, à medida que a energia aumenta, os valores interpolados de t tendem a oscilar em torno dos resultados experimentais, pois os efeitos de absorção fotoelétrica para altas energias tornam-se desprezíveis em relação às outras formas de absorção da energia do fóton, e têm comportamento de díficil descrição. Por fim, a seção de choque total, integrada em um determinado grupo g é

(14)

Tabela 1
- Valores numéricos das seções de choque microscópicas de absorção fotoelétrica (t), espalhamento Compton (s) e produção de pares (k), obtidas dos dados de Davisson & Evans (1952) (DE), de regressão linear múltipla (RL) e da distribuição de Klein-Nishina (KN).

ENERGIA (Kev)	t(DE) (barns)	s(DE) (barns)	k(DE) (barns)	t(RL) (barns)	s(KN) (barns)	k(RL) (barns)
102,2	0,057	0,490	_	0,0572	0,4716	_
340,1	0,00127	0,335	_	0,00122	0,3327	_
681,1	0,000184	0,253	_	0,000177	0,2485	_
1022,0	< 10-4	0,209	_	< 10-4	0,2057	_
2043,0	< 10-5	0,145	0,002225	< 10-5	0,1411	0,00227
4086,0	<10-6	0,095	0,01008	< 10-6	0,0949	0,00993

Table 1 - Photoelectric absorption (t), Compton scattering (s) and pair production (k) numerical microscopic cross-section values obtained from the Davisson & Evans data (DE), multiple linear regression scheme (RL) and the Klein-Nishina distribution (KN).

SIMULAÇÃO DE FONTES ALEATÓRIAS

As principais fontes de radiação gama natural da Terra são devidas, quase exclusivamente, aos radioisótopos K-40, U-238 e Th-232 (Meyers, 1992). Vários autores (Serra, 1984; Willie, 1984; entre outros) determinaram através de experimentos de laboratório os teores destes radioelementos, estimando suas concentrações mínimas e máximas (C_K,min e C_K,max, C_U,min e C_U,max, C_Th,min e C_Th,max) nos principais minerais formadores de rochas (Wahl,1983). Os resultados deste trabalho estão mostrados na Tab. 2.

Tabela 2
- Concentrações mínimas e máximas de K-40, U-238 e Th-232 nos minerais.

MINERAL	CK,min- CK,max	CU,min- CU,max	CTh,min- CTh,max
Quartzo	0,00 - 0,15	0,00 - 0,40	0,00 - 0,20
Calcita	0,00 - 0,10	0,00 - 1,00	0,00 - 0,50
Dolomita	0,00 - 0,40	0,00 - 1,00	0,00 - 1,50
Anidrita	0,00 - 0,10	0,00 - 0,35	0,00 - 0,10
Gipsita	0,00 - 0,15	0,00 - 1,20	0,00 - 1,00
Caolinita	0,00 - 0,42	1,50 - 3,00	6,00 - 19,00
Ilita	0,00 - 4,50	0,00 - 1,50	6,00 - 19,00
Montmorilonita	0,00 - 0,16	2,00 - 5,00	14,00 - 2,.0
K-Feldspato	11,80 - 14,00	0,00 - 1,50	0,00 - 0,10
Muscovita	7,90 - 9,80	0,00 - 1,00	0,00 - 1,00
Biotita	6,70 - 8,30	0,00 - 1,00	0,00 - 1,00
Água	0,00 - 0,00	0,00 - 0,00	0,00 - 0,00

Table 2 - Minimum and maximum concentrations of K-40, U-238 and Th-232 in minerals.

O objetivo desta seção é determinar as concentrações em peso desses radioelementos supondo-os distribuídos de forma aleatória ao longo das camadas, de acordo com certas funções de distribuição de probabilidade (fdp) e tendo por base os valores mínimos e máximos obtidos experimentalmente (Tab. 2). Para tal, faz-se uso do método Monte Carlo para gerar aleatoriamente as fontes gama naturais dentro de cada elemento da malha de elementos finitos. Obviamente, estas fontes estarão dispostas de acordo com as diferentes regiões dos modelos. A fdp terá, assim, caráter predominantemente geológico e apresentará valores médios , e , das concentrações de potássio, urânio e tório, distribuídas uniformemente da forma a seguir:

, (15)

, (16)

, (17)

onde r1, r2 e r3 são números pseudo-aleatórios distribuídos uniformemente no intervalo [0,1]. Os conteúdos efetivos aleatórios de cada isótopo, C_K, C_U e C_Th, serão simulados através da distribuição gaussiana, tendo como médias as concentrações das Eqs. (15), (16) e (17) e desvios padrão , e .Neste ítem, será adotado o método de rejeição descrito por Rubinstein (1981) e Kalos & Whitlock (1986). Se para cada concentração dentro de um determinado elemento da malha FEM for observado que C_K,min£ C_K£ C_K,max, C_U,min£ C_U £ C_U,max e C_Th,min£ C_Th£ C_Th,max , as três concentrações serão aceitas. Caso alguma das condições impostas não seja observada, a concentração a ela relacionada será rejeitada e sua simulação refeita quantas vezes forem necessárias, e sempre a partir de um novo número pseudo-aleatório, até que a condição previamente imposta seja definitivamente observada.

DISCUSSÃO

Para testar a eficiência do método, os fluxos escalares obtidos com o algoritmo FEM foram comparados com os espectros calculados por Bertozzi et al. (1981), considerando modelos com fontes uniformemente distribuídas. As soluções de Bertozzi et al. (1981) consideram que o espectro é o resultado de um fenômeno local, originado pela degradação dos fótons de alta energia devida ao espalhamento Compton, à absorção fotoelétrica e ao processo de formação de pares no ponto de interesse. A seguir, um teste adicional foi feito comparando a resposta da modelagem FEM para G = 60 grupos de energia com respostas para G = 10, 20, 30, ... grupos, sabendo-se de antemão que 60 grupos eram suficientes para atingir a convergência do fluxo escalar. Os resultados mostraram que para G > 46 grupos, nenhum ganho adicional em termos de convergência pode ser conseguido e este é o valor de G adotado para as simulações desenvolvidas neste trabalho.

Modelo de calibração

Uma etapa importante deste trabalho consiste em estabelecer uma relação direta entre a taxa de contagem total, ou seja, o fluxo angular, e a distribuição e quantidade de material radioativo na formação. Isto é feito através da calibração da sonda em um poço com o mesmo padrão das companhias de petróleo, isto é, 20 cm de diâmetro, preenchido com água de densidade 1000 kg/m³, porém com a sonda centralizada para atender a geometria adotada neste trabalho. A Fig. 4 mostra os espectros total e parciais dos três isótopos considerados e as janelas utilizadas para suas detecções. A janela do potássio está centrada na energia 1460,0 keV e corresponde ao radioisótopo K-40 (Figs. 4(a) e 4(b)). A janela do urânio está centrada na energia 1765,0 keV e corresponde ao pico do Bi-214 pertencente à série radioativa natural do U-238 (Figs. 4(a) e 4(c)). A janela do tório está centrada na energia 2615,0 keV e corresponde ao pico do Tl-208 pertencente à série do Th-232 (Figs. 4(a) e 4(d)). O modelo hipotético que originou o espectro total é um arenito com concentrações C_K = 4,0%, C_U = 12,0 ppm e C_Th = 24,0 ppm e que reproduz a formação artificialmente radioativa, construída na Universidade de Houston (Ellis, 1987). Os detalhes destes modelos de calibração são descritos em Belknap et al. (1976). Para os espectros parciais, os modelos são litologicamente semelhantes, mas as concentrações em massa dos radioisótopos são (1) espectro do potássio: C_K = 4,0%, C_U = C_Th = 0; (2) espectro do urânio: C_U = 12,0 ppm, C_K = C_Th = 0, (3) espectro do tório: C_Th = 24,0 ppm C_K = C_U = 0. É a calibração em formações padrões que permite determinar as concentrações em massa destes radioisótopos.

Figura 4
- Perfis sintéticos para calibração da sonda: (a) espectro total, (b) espectro do potássio-40, (c) espectro do urânio-238 e (d) espectro do tório-232.

Figure 4 - Synthetic logs for sonde calibration: (a) total spectrum, (b) potassium-40 spectrum, (c) uranium-238 spectrum and (d) thorium-232 spectrum.

O método de análise desenvolvido aqui tem como base a combinação dos modelos propostos por Koizumi (1988) e Smith et al. (1983), para calibração e correções ambientais de um sistema geral de perfilagem usando raios gama naturais. Como as concentrações em massa [C_K, C_U , C_Th] dos radioelementos presentes na formação são estipuladas a partir das taxas de contagem [r_K, r_U , r_Th] em cada janela, e que estas relações são lineares, tem-se o seguinte sistema:

(18)

onde cada elemento a_ij da matriz de sensibilidade é uma constante determinada diretamente dos modelos de calibração, para um arranjo geométrico específico da sonda. Cada elemento a_ij representa a contribuição do radioelemento j na i-ésima janela de detecção. Como exemplo, citamos o caso genérico de um modelo com as seguintes concentrações : C_K = 4,0%, C_U = C_Th = 0. Levando estes valores à Eq. (18), as constantes referentes às contribuições da concentração de potássio nas três janelas selecionadas serão dadas pelas relações a₁₁= r_k/4,0, a₂₁= r_u/4,0 e a₃₁= r_th/4,0 . No entanto, a matriz de sensibilidade A pode ser obtida completamente, perfilando M ³ 3 modelos de calibração a fim de se obter uma matriz resposta R, cujos elementos da j-ésima coluna são as taxas de contagem devidas ao modelo j, nas três janelas de detecção,

(19)

ou de uma forma mais compacta,

R=AC. (20)

A solução deste sistema será obtida pós-multiplicando ambos os lados da Eq. (20) por C^T,

RC^T = ACC^T (21)

a partir da qual se pode determinar a matriz de sensibilidade invertendo o sistema anterior,

A = (RC^T)(CC^T)^-1 . (22)

De posse da matriz de sensibilidade, as concentrações C_K , C_U e C_Th , para quaisquer modelos, podem ser determinadas pela seguinte expressão

c = T r, (23)

onde a matriz de transformação T é dada por

T = (CC^T)(RC^T)^-1 . (24)

No entanto, sabe-se que os resultados obtidos usando a Eq. (23) são severamente afetados por flutuações estatísticas, com exceção do tório. Na prática, estes efeitos são minimizados por diversas técnicas de filtragem, sendo a principal delas, a que utiliza o filtro de Kalman. Sabe-se ainda que grande parte dos erros de calibração e interpretação ocorrem porque as concentrações nos modelos de calibração são medidas em um único ponto dentro do poço, diferente dos trabalhos correntes de perfilagem onde a sonda é deslocada continuamente ao longo do poço enquanto se fazem as medições. Ao assumir uma relação linear entre as taxas de contagens e as concentrações em massa dos radioisótopos mencionados, supõe-se que os elementos de interesse estão uniformemente distribuídos em volta da sonda e em concentrações relativamente altas. Isso porém não ocorre na prática, pois aí se verificam concentrações bastante reduzidas e distribuídas de forma anômala. Além disso, as respostas são consideradas não-lineares por causa da degradação da energia dos raios gama nas janelas de menor energia. Neste caso, é adotada uma técnica de correção conhecida como "stripping". Como a concentração de tório é a única a ser determinada satisfatoriamente, devido à independência de sua janela em relação aos outros radioisótopos, C_Th será mantido e utilizado para corrigir as outras concentrações. Isto quer dizer que se deve manter inalterado o elemento t₃₃ da matriz de transformação T.

Seja f (x, y) um fator que determina a contribuição exclusiva da concentração do radioelemento y na janela do radioelemento x. Sejam r_(K,Th), r_(U,Th) e r_(Th,Th), as taxas de contagem nas janelas de potássio, urânio e tório, devidas à concentração de tório C_Th = 24,0 ppm. Os fatores de contribuição f (K, Th) e f (U, Th) são determinados da seguinte forma:

cps/ppm de tório, (25)

cps/ppm de tório. (26)

A seguir repete-se a tarefa anterior para medir as taxas de contagem r_(K,U) e r_(U,U), devidas à concentração de urânio (C_U = 12,0 ppm). Assim, pode-se determinar da mesma forma o fator f (K, U) da contribuição do urânio na janela de potássio:

cps/ppm de urânio. (27)

A partir destes fatores de contribiução pode-se finalmente corrigir as taxas de contagem r_K e r_U medidas pela sonda de perfilagem. As taxas corrigidas serão denotadas por r’_K e r’_U e assumirão os seguintes valores:

(28)

(29)

De posse das taxas corrigidas, pode-se calcular as concentrações em massa definitivas de potássio (C_K) e urânio (C_U). Para isto, mede-se r(K,K), a taxa de contagem na janela de potássio devido a uma concentração em massa C_K = 4,0. Os fatores t_K e t_U para a transformação das taxas r’_K e r’_U em concentrações de potássio e urânio, são

(30)

(31)

Estas considerações permitem reconstruir a matriz T para as taxas de contagem corrigidas, supondo que , ou seja, o urânio apresenta contribuição desprezível na janela do tório:

(32)

É importante lembrar que os valores numéricos obtidos com as análises acima, só se aplicam ao modelo idealizado neste trabalho. Outros casos devem ser revistos a fim de produzirem resultados condizentes com os modelos geológicos propostos.

Modelos sintéticos

Na Fig. 5, são mostrados os perfis das concentrações em massa de (a) potássio-40, (b) urânio-238 e (c) tório-232 para um modelo constituído de nove camadas plano-paralelas com as especificações descritas na Tab. 3. Um primeiro detalhe observado nestes perfis é a reversão nas amplitudes de tório e urânio com relação ao potássio para as camadas c2 e c3. Este detalhe, mesmo em face dos pequenos contrastes entre as concentrações dos três radioisótopos, permite identificar as profundidades das interfaces que limitam as camadas c1, c2 e c3, o que certamente não seria possível sem esta reversão. Nota-se que, mesmo apresentando contrastes apreciáveis, as camadas c6, c7 e c8 são completamente indistinguíveis. Neste caso, a indistinguibilidade deve-se à reduzida espessura destas camadas, aquém da resolução da sonda, cujo valor é de 30 cm e corresponde à altura do cristal NaI(Tl), de forma cilíndrica e posicionado verticalmente no interior da sonda para permitir incidência lateral do fluxo de radiação gama natural (Ellis, 1987). A fim de proceder-se a uma análise predominantemente numérica dos resultados das Figs. 5(a), (b) e (c), tomemos como exemplo três pontos centrais da camada c2, localizados nas profundidades de 80, 85 e 90 cm. Os valores mínimo e máximo possíveis para as concentrações desta camada, baseado nas Tabs. 2 e 3, são C_K,min = 0% e C_K,max = 1,413%; C_U,min = 0,6ppm e C_U,max = 1,99ppm; C_Th,min = 6,0ppm e C_Th,max = 12,92ppm, com as correspondentes médias aritméticas C_K(80,85,90) = 0,7065%, C_U(80,85,90) = 1,295ppm e C_Th(80,85,90) = 9,46ppm. A aleatoriedade do método utilizado resultou em perfis com diferentes valores de C_i,j simulados pelo método FEM, onde i representa o radioisótopo e j a profundidade. Para os pontos em análise temos: C_K,80 = 0,5583%, C_K,85 = 0,5702%, C_K,90 = 0,5877% e média C_K,M = 0,5721%; C_U,80 = 1,2229ppm, C_U,85 = 1,2198ppm, C_U,90 = 1,2187ppm e média C_U,M = 1,2205ppm; C_Th,80 = 9,1689ppm, C_Th,85 = 9,2501ppm, C_Th,90 = 9,2853ppm e média C_Th,M = 9,2348ppm. Vê-se que estes valores médios simulados diferem em 23,5%, 6,1% e 2,45% das médias aritméticas dos valores mínimo e máximo possíveis para a camada, para as concentrações de potássio, urânio e tório, respectivamente. Embora dentro das variações estipuladas pelo modelo estatístico empregado (distribuição de Poisson), as maiores discrepâncias estão relacionadas com as camadas de menor espessura, o que implica em deflexões mais acentuadas nos perfis já a partir de seus pontos mais centrais. Esta mesma análise quando refeita para a camada espessa c5, nas profundidades de 320, 325 e 330 cm, mostrou variações de 0,18%, 9% e 0,39% para as concentrações de potássio, urânio e tório. É importante ressaltar que de uma forma geral, os processos de "stripping" têm efeito predominante sobre estes resultados. Por exemplo, a menor discrepância para uma mesma camada é observada para as concentrações de tório que são praticamente isentas destas correções. Por outro lado, vê-se que o urânio é realmente o menos confiável devido em parte ao seu complexo comportamento geoquímico, e em parte à localização de sua janela de detecção.

Tabela 3
- Especificações gerais do modelo da Fig. 5: camadas (CAM.), espessura (ESP.), profundidades do topo e base (TOPO - BASE) e frações mineralógicas (FRAC.MIN.). Nomenclatura: QTZ(quartzo), CAL(calcita), DOL(dolomita), ANI(anidrita), GIP(gipsita), CAO(caolinita), ILI(ilita), MON(montmorilonita), FEL(K-feldspato), BIO(biotita), MUS(muscovita), f (porosidade com SW=100 %).

CAM.	ESP. (cm)	TOPO-BASE (cm-cm)	FRAC.MIN.
c1	50	0-50	QTZ=50% ; GIP=20% ; ILI=10% ; f=20% ;
c2	70	50-120	QTZ=10% ILI=30% ; MON=30% ; f=30% ;
c3	60	120-180	DOL=20% BIO=20% ; CAO=30% ; f=30% ;
c4	50	180-230	QTZ=50% MON=10% ; ANI=20% ; f=20% ;
c5	90	230-420	QTZ=10% FEL=30% ; MON=30% ; f=30% ;
c6	30	420-450	QTZ=20% MON=20% ; CAO=30% ; f=30% ;
c7	40	450-490	QTZ=50% ILI=10% ; CAL=20% ; f=20% ;
c8	30	490-520	QTZ=10% DOL=30% ; CAL=30% ; f=30% ;
c9	80	520-600	QTZ=20% MUS=20% ; MON=30% ; f=30% ;
POÇO	600	0-600	SW=100%

Table 3 - General specifications for the model of Fig. 5: Layer (CAM.), thickness (ESP.), top and bottom depth (TOPO - BASE) and mineralogical fractions (FRAC.MIN.). Nomenclature: QTZ(quartz), CAL(calcite), DOL(dolomite), ANI(anhydrite), GIP(gypsite), CAO(kaolinite), ILI(illite), MON(montmorillonite), FEL(K-feldspar), BIO(biotite), MUS(muscovite), f (porosity with SW=100 %).

Figura 5
- Perfis sintéticos de espectrometria de raios gama naturais, mostrando as concentrações mínima (linha tracejada), máxima (linha cheia) e simuladas (linha cheia-pontuada) de cada camada: (a) concentração em massa de pótassio-40, (b) concentração em massa de urânio-238, (c) concentração em massa de tório-232.

Figure 5 - Natural gamma-ray spectrometry synthetic logs showing the minimum (dashed line), maximum (solid line) and simulated (solid-pointed line) concentrations by weight of (a) potassium-40, (b) uranium-238 and (c) thorium-232.

Na Fig. 6, são também mostrados os perfis das concentrações em massa de (a) pótassio-40, (b) urânio-238 e (c) tório-232 para um modelo constituído de 19 camadas plano-paralelas com as especificações descritas na Tab. 4. Na camada c4, os baixos níveis de radiação se devem à ausência de minerais de argila, feldspatos e micas que detêm a maior parte da radioatividade natural dos sedimentos (Tab. 1). As camadas c6, c7 e c8 são indeterminadas se se considerar apenas o perfil de urânio (Fig 6(b)). Isto dá inclusive a falsa idéia de existir entre as profundidades de 350,0 e 520,0 cm uma só camada espessa. Isso é, no entanto, desmentido ao se observar os perfis de potássio e urânio, embora ainda restem dúvidas quanto à existência da camada c8, pois os baixos contrastes nas concentrações das camadas c7 e c8 não permitem visualisar direito a existência da interface entre elas. As camadas c12 e c13 por serem estremamente delgadas não são distinguidas através do perfil de potássio. No entanto, devido ao contraste considerável nas concentrações de urânio e tório, pode-se distinguir aí a ocorrência de uma única camada englobando c12 e c13. Outras duas camadas adjacentes que ocorrem indistinguíveis neste modelo são as camadas c16 e c17, embora possam ser identificadas como uma única estrutura nos perfis das Figs. 6(a), 6(b) e 6(c). Uma observação importante neste modelo ocorre no perfil de tório (Fig. 6(c)). Embora este perfil seja reconhecidamente o de melhor resolução de litologias, nota-se que nas camadas c18 e c19 não estão evidenciadas duas camadas, mas apenas uma, a despeito do alto contraste entre elas. Esta observação no entanto tem uma explicação simples. Como os valores mínimos e máximos das concentrações de tório nestas duas camadas perfazem a mesma média, obviamente os valores de CTh para ambas serão em média iguais. Repetindo-se para o modelo da Fig. 6 o mesmo critério de análise numérica do modelo da Fig. 5, obteve-se discrepâncias entre as médias aritméticas dos valores mínimo e máximo possíveis da camada c1 e seus respectivos valores médios simulados via método FEM de 20,5% para o potássio, 17,8% para o urânio e 3,95% para o tório, corroborando a maior confiabilidade na medição da concentração deste último no ponto de estudo.

Tabela 4
- Especificações gerais do modelo da Fig. 6: camadas (CAM.), espessura (ESP.), profundidades do topo e base (TOPO - BASE) e frações mineralógicas (FRAC.MIN.). Nomenclatura: QTZ(quartzo), CAL(calcita), DOL(dolomita), ANI(anidrita), GIP(gipsita), CAO(caolinita), ILI(ilita), MON(montmorilonita), FEL(K-feldspato), BIO(biotita), MUS(muscovita), f (porosidade com SW=100 %).

CAM.	ESP. (cm)	TOPO-BASE (cm-cm)	FRAC.MIN.
c1	80	0-80	ANI=50% ; GIP=20% ; ILI=10% ; f=20% ;
c2	50	80-130	DOL=10% ; CAO=30% ; ILI=30% ; f=30%
c3	50	130-180	CAL=20% ; CAO=30% ; BIO=20% ; f=30%
c4	120	180-300	QTZ=50% ; CAL=20% ; DOL=10% ; f=20%
c5	50	300-350	QTZ=10% ; FEL=30% ; MON=30% ; f=30%
c6	100	350-450	QTZ=20% ; FEL=20% ; ILI=30% ; f=30%
c7	40	450-490	QTZ=50% ; ILI=10% ; CAL=20% ; f=20%
c8	304	90-520	QTZ=10% ; DOL=30% ; CAL=30% ; f=30%
c9	80	520-600	QTZ=20% ; CAO=30% ; MUS=20% ; f=30%
c10	50	600-650	QTZ=50% ; ANI=20% ; GIP=10% ; f=20%
c11	70	650-720	QTZ=10%; MON=30%; ILI=30%; f=30%
c12	10	720-730	DOL=20%; CAO=30%; BIO=20%; f=30%
c13	10	730-740	QTZ=50%; ANI=20%; FEL=10%; f=20%
c14	60	740-800	QTZ=10%; CAO=30%; FEL=30%; f=30%
c15	50	800-850	QTZ=20%; CAO=30%; MUS=20%; f=30%
c16	40	850-890	QTZ=50%; GIP=20%; ILI=10%; f=20%
c17	30	890-920	QTZ=10%; DOL=30%; CAL=30%; f=30%
c18	30	920-950	QTZ=20%; MON=30%; MUS=20%; f=30%
c19	50	950-1000	QTZ=20%; CAO=30%; ILI=20%; f=30%
POÇO	1000	0-1000	SW=100%

Table 4 - General specifications for the model of Fig. 6: Layer (CAM.), thickness (ESP.), top and bottom depth (TOPO - BASE) and mineralogical fractions (FRAC.MIN.). Nomenclature: QTZ(quartz), CAL(calcite), DOL(dolomite), ANI(anhydrite), GIP(gypsite), CAO(kaolinite), ILI(illite), MON(montmorillonite), FEL(K-feldspar), BIO(biotite), MUS(muscovite), f (porosity with SW=100 %).

Figura 6
- Perfis sintéticos de espectrometria de raios gama naturais, mostrando as concentrações mínima (linha tracejada), máxima (linha cheia) e simuladas (linha cheia-pontuada) de cada camada: (a) concentração em massa de pótassio-40, (b) concentração em massa de urânio-238, (c) concentração em massa de tório-232.

Figure 6 - Natural gamma-ray spectrometry synthetic logs showing the minimum (dashed line), maximum (solid line) and simulated (solid-pointed line) concentrations by weight of (a) potassium-40, (b) uranium-238 and (c) thorium-232.

CONCLUSÕES

Um modelo teórico baseado nos métodos de elementos finitos e Monte Carlo foi desenvolvido para simular perfis sintéticos de espectrometria de raios gama naturais em meios aleatórios. Ao desprezar a aplicação direta da lei de Fick como ferramenta tradicional na redução da equação de transporte de Boltzmann, introduzindo em seu lugar, a expansão em harmônicos esféricos do fluxo angular e da fonte, o fluxo escalar na região de altas energias pôde ser obtido de forma mais eficaz. Verificou-se que os modelos do sistema K-U-Th são, na verdade, sensíveis à ação conjunta da quantidade e da distribuição do material radioativo nas rochas. Nos modelos hipotéticos de calibração, foram observadas fontes distribuídas uniformemente e simulando de duas a cinco vezes mais a radioatividade de um folhelho típico, o que implica numa maior profundidade de investigação, capaz de não refletir exatamente as situações reais. Contudo, outros modelos de calibração com menores índices de radioatividade foram propostos e combinados com o modelo padrão a fim de se conseguir níveis de estatísticas aceitáveis para os modelos apresentados. Os perfis sintéticos mostraram que análises baseadas apenas nas concentrações em massa de pótassio-40 e o tório-232 definem de forma mais eficaz as interfaces de transição entre diferentes litologias, desde que suas espessuras sejam maiores do que 30 cm, que é a resolução da sonda simulada neste trabalho. O urânio, por apresentar maior capacidade de migração e carreamento pelos fluidos que percolam as rochas, devido à suas características geoquímicas já tão bem conhecidas e amplamente divulgadas, não se constitui em um indicador eficaz de litologia, fato este já evidenciado por diversos trabalhos geoquímicos. Enfim, fica evidenciado que os baixos contrastes de concentração entre camadas adjacentes, principalmente para o urânio, seguido do potássio e em menor escala o tório, não permitem determinar satisfatoriamente as intricadas subdivisões de litologias complexas.

AGRADECIMENTOS

O autor agradece o apoio financeiro recebido através do CNPq, Ref.: Proc.: 30.0473/95.8 (NV) - Desenv. Cient. Regional, e às facilidades computacionais do Departamento de Geofísica da UFPa, do NCE-Núcleo de Computação Eletrônica da UFRJ e do Setor de Geologia de Engenharia e Ambiental do Departamento de Geologia da UFRJ, que tornaram este trabalho possível.

Received: May, 1996

Accepted: November, 1998

SIMULATION OF NATURAL GAMMA-RAY SPECTOMETRY LOGS IN RANDOM MEDIA

Synthetic natural gamma-ray spectrometry logs are simulated for radioisotopes sources K-40, U-238 and Th-232 randomly distributed in formations. In this case, the transport of gamma radiation is calculated by the finite-element method (FEM) as 2 1/2-dimensional scalar flux from the multigroup diffusion approximation of the Boltzmann transport equation in phase space (Nowosad, 1978). This solution is known as Pl approximation, where the variables direction and random spectral sources are expanded in spherical harmonics through the orthogonal Legendre polynomials (Etherington, 1958). The isotopic sources are simulated by a simple process of random source generation based on the Monte Carlo method (Rubinstein, 1981). In this method the physical parameters of formation are established by means of nuclear and non-relativistic physics concepts according to Klein-Nishina distributions (Davisson & Evans, 1952). The concentrations by weight of K-40, U-238 and Th-232 are not constant for each lithology. They assume values varying randomly within each element of the FEM mesh concerning the variation interval of each radioisotope for a specific mineral (Wahl, 1983). Next, the proposed methodology is applied to obtain synthetic well logs and to show the multiple applications of this important logging system in complex lithologies. To apply accurately these ideas, there must be a valid geological and mathematical model of the physics of the system under study. The idealized geological model is shown in Fig. 1, while the ideal mathematical model is that suggested by the Boltzmann transport equation. Hence to get a solution by the FEM, Eq. (1) must be reduced by the P₁ approximation. This reduction is performed having in mind that the flux depends only on the spatial coordinate and physical parameters of the formation. An accurate system based on models proposed by Koizumi (1988) and Smith et al. (1983) are developed for the calibration of the hypothetic sonde to be simulated in this paper. The spectral response of the sonde are shown in Fig. 4. This calibration model allows us to determine CK, CU and CTh by means of Eq. (23) through the transformation matrix given in Eq. (32). In Fig. 5 an important inversion on the thorium and uranium concentration values is observed in layers c2 and c3 related to that of potassium. This detail allows us to define accurately the boundary interfaces of the layers described in Tab. 3, which could be impossible without this inversion of data. It can be seen that despite the high contrast between their concentrations, the layers c6, c7 and c8 are undistinguishable mainly because the thicknesses of these thin layers are bellow the sonde resolution (30 cm). In Fig 6, a model constituted by nineteen layers is shown. The specifications of this model are given in Tab. 4. The low radiation observed in layer c4 are due to the absence of clays, feldspars and mica minerals. Layers c6, c7 and c8 are not distinguishable if only the uranium log is considered (Fig. 6(b)) and give us the impression of a thick layer extending from 350.0 to 520.0 cm observing the potassium and thorium logs (Figs. 6(a) and 6(c)).

NOTES ABOUT THE AUTHOR

NOTAS SOBRE O AUTOR

O autor é graduado em Engenharia Civil e de Transportes, tendo se especializado em Ciências e Técnicas Nucleares pelo Departamento de Engenharia Nuclear da UFMG. Concluiu o doutorado em geofísica em 1993, com ênfase em Geofísica de Poço, pelo Departamento de Geofísica da UFPa. Exerce atualmente o cargo de Professor Adjunto no Setor de Geologia de Engenharia e Ambiental do Departamento de Geologia da UFRJ. Sua área de interesse é a simulação numérica de perfis nucleares de poços através dos métodos de Elementos Finitos e Monte Carlo.

• BARRY, J. M. & POLLARD, J. P.-1982- Solution of neutron diffusion linear equations by the method of implicit non-stationary iteraction. North-Holland Publishing, New York.
• BELKNAP, W. B., DEWAN, J. T., KIRKPATRICK, C. V., MOTT, W. E., PEARSON, A. J. & RABSON, W. R.-1976- API calibration facility for nuclear logs. Drill Prod. Prac., API, Houston.
• BERTOZZI, W., ELLIS, V. & WAHL, J. S.-1981- The physical foundation of formation lithology logging with gamma rays. Geophysics, 46: 120 - 132.
• DAVIES, A. L.-1980- The finite element method: a first approach. Claredon Press. 287pp.
• DAVISSON, C. M. &, EVANS, R. D.-1952- Gamma-ray absorption coefficients. Review of Modern Physics, 24: 79-107.
• ELLIS, D. V.-1987- Well logging for earth scientists. New York, Elsevier, 532pp.
• ETHERINGTON, H.-1958- Nuclear Engineering Handbook, McGraw-Hill, New York. 1012pp.
• KALOS, M. H. & WHITLOCK, P. A.-1986- Monte Carlo Methods. John Wiley & Sons, New York. 186pp.
• KOIZUMI, C. J.-1988- Computer determination of calibration and environmental corrections for a natural spectral gamma ray logging system. SPE Formation Evaluation, 20: 637 - 644.
• MEYERS, G. D.-1992- A review of nuclear logging. The Log Analyst, 34: 228-238.
• NOWOSAD, P.-1978- Elementos da teoria dos reatores nucleares, In: Nowosad, P., Operadores positivos e optimizaçăo - aplicaçăo á engenharia nuclear. CBPF-Centro Brasileiro de Pesquisas Físicas, Rio de Janeiro.
• ODEN, J. T., & REDDY, J. N.-1976- An introduction to the mathematical theory of finite elements. John Wiley & Sons. 429pp.
• RUBINSTEIN, R. Y.-1981- Simulation and Monte Carlo Method. New York, John Wiley & Sons, 278pp.
• SERRA, O.-1984- Fundamentals of well-log interpretation. Amsterdam, Elsevier, 423pp.
• SMITH, H. D., REBBINS, C. A., GARDNER, L. L. & DEATON, J. G. -1983- A multifunction compensated spectral natural gamma ray logging system. In: Annual technical conference and exhibition, San Francisco, Proceeding: SPE.
• WAHL, J. S.-1983- Gamma-ray logging. Geophysics, 48: 1536-1550.
• WILLIE, A. W.-1984- Nuclear assaying of mining borehole: an introduction. Amsterdam, Elsevier, 344pp.

Datas de Publicação

• Publicação nesta coleção
  23 Dez 1999
• Data do Fascículo
  Jul 1998

Histórico

• Recebido
  Maio 1996
• Aceito
  Nov 1998