基于多尺度震源模型的1905年康格拉地震确定性宽频带地震动模拟研究

本研究针对1905年Mw7.8康格拉地震缺乏仪器记录数据的问题，通过开发西喜马拉雅地区三维谱元计算模型，采用多尺度震源模型成功模拟了0.001-20 Hz宽频带地震动。研究验证了模型在0.5 Hz以下的准确性，并通过与宏观烈度记录的对比证实了模拟结果的可靠性。该研究为喜马拉雅地区地震危险性评估提供了重要技术框架，对工程结构抗震设计具有重要指导意义。

　　

1905年4月4日发生的M_w7.8康格拉地震是喜马拉雅地震史上最具破坏性的地震之一，造成了广泛破坏，麦卡利烈度达到X级（极端）。由于当时缺乏强震仪器记录，这个重大历史地震的许多关键特征至今仍未被充分认识。地震导致包括旁遮普、喜马偕尔邦和查谟克什米尔在内的大面积区域遭受严重破坏，约10万座建筑被毁，2万人丧生。除了人员伤亡，地震还严重损坏了许多历史建筑，如康格拉堡垒和马苏鲁岩石切割寺庙，并引发了大规模山体滑坡。

由于缺乏当时的地震动记录，对该地震的认识主要来自印度地质调查局（GSI）的第一手观测资料。早期记录最初估计地震震级为M_w≥8，等震线图显示距震中250公里的德拉敦的罗西-福雷尔（RF）烈度为VIII度。然而，多项研究重新评估了该事件的规模和宏观烈度，发现初始估计被夸大了。后来，该事件的震级被确定为M_w7.8，烈度图按梅德韦杰夫-斯蓬霍伊尔-卡尼克（MSK）烈度标度进行解释。

康格拉地震最引人关注的方面之一是地表破裂，由于滑坡特别是在等震区的山区，极难绘制。关于该事件的破裂特征也存在不确定性，地震被认为起源于主中央逆冲断层并延伸至喜马拉雅前缘逆冲断层。有研究推断破裂沿贾瓦拉穆克希逆冲断层发生，而另一些研究则利用卫星图像和野外填图表明事件沿康格拉谷断层（KVF）破裂，这是位于主边界逆冲断层和贾瓦拉穆克希逆冲断层之间的一条活动断层。

自1905年康格拉事件以来，喜马拉雅地区发生了许多重大地震（M_w>7），如1934年比哈尔-尼泊尔地震（M_w8.1）、1947年阿鲁纳恰尔邦地震（M_w7.9）、1950年阿萨姆地震（M_w8.6）、2005年穆扎法尔布尔地震（M_w7.6）和2015年尼泊尔地震（M_w7.9）。然而，除2015年尼泊尔地震外，这些事件均未被记录，原因相同。因此，多项研究使用多种方法模拟了印度历史地震的地震动。

一些研究采用经验格林函数（EGF）技术模拟高破坏性历史地震，如2001年库奇地震和1905年康格拉地震。与使用经验格林函数（EGF）技术模拟1905年康格拉地震的研究不同，本研究采用完全基于物理的确定性框架，将高分辨率三维谱元模型与多尺度震源表征相结合。虽然EGF方法依赖于较小地震记录的线性缩放，通常限制频率范围并忽略复杂的源-路径相互作用，但本方法产生连续的宽频带地震动（0.001-20 Hz），无需频率拼接，明确包含三维速度非均质性、地形和盆地效应。这能够真实再现低频和高频分量，改进PGA分布的空间细节（特别是在恒河盆地），并与宏观烈度观测更接近。这些能力增强了模拟的物理真实性和工程相关性，解决了早期基于EGF研究的关键局限性。

此外，多项研究采用随机有限断层地震学模型生成强震动时程，如1897年西隆地震、2005年克什米尔地震和1950年阿萨姆地震，并结合分析方法。随机地震学模型在模拟高频地震动（约1 Hz以上）方面最有效，但难以准确模拟低频运动（<1 Hz），这对评估大型结构和长周期地震动效应至关重要。此外，场地特定效应使用放大系数或通用场地响应函数建模，可能无法准确表示当地场地条件。事实上，EGF和随机有限断层方法都不能准确考虑材料特性的变化、喜马拉雅地形和印度恒河平原。因此，使用有限元方法进行基于物理的建模，开发高自由度计算模型并生成地震动。然而，基于确定性方法生成的地震动频率限制在1至2 Hz。因此，多项研究结合确定性和随机方法来生成宽频带地震动。最近，有研究通过结合使用谱元法获得的低频波形和使用改进的随机地震学模型生成的高频地震图，模拟了1934年比哈尔-尼泊尔地震的宽频带地震动。然而，结合使用不同方法生成的地震动对目标频率和带宽的选择高度敏感。

随着超级计算能力的全球革命，许多研究开始使用纯基于物理的确定性数值模型生成主要地震的宽频带地震动。因此，多项研究使用多尺度或复合震源模型来生成宽频带地震动。此外，随机滑动分布广泛用于地震震源建模以捕捉破裂变异性，这与模拟空间和参数变异性的不确定性量化进展并行。最近的工作表明，将空间变异性建模与深度学习代理模型耦合可提高大规模模拟的计算效率。然而，许多随机方法在完全纳入断层带复杂性和保持计算可行性方面面临限制，推动了在宽频带地震动模拟中平衡物理真实性和效率的框架的使用。有研究将多尺度震源模型实施到谱元法中，模拟1994年北岭地震的地震动。在喜马拉雅事件的背景下，迄今为止确定性数值模型仅包括简单的震源描述，即小地震的点源和大地震的有限断层模型。然而，简单模型不足以精确模拟高频（对低层结构至关重要）和低频（对大型结构重要）地震动。此外，近断层区域观察到的永久地面位移、上盘效应和抛掷步无法用简单的断层表征充分建模。除了震源特征外，模拟具有材料特性变化的介质以及IGP和喜马拉雅地形的影响具有挑战性，特别是在没有记录数据验证模拟地震动的情况下。

因此，本研究重点开发西喜马拉雅地区的三维计算模型，并使用2013年M_w5.8巴达瓦地震记录的速度时程验证模拟地震动。经验证后，可以模拟1905年康格拉地震的低频速度时程。此外，研究旨在开发1905年康格拉地震的多尺度震源模型，以确定性模拟近断层区域的宽频带地震动。在缺乏该历史事件记录地震动的情况下，可以将报告的宏观烈度值与模拟地震动的烈度值进行比较。

研究人员主要采用了谱元法（Spectral Element Method, SEM）进行三维地震动模拟，开发了西喜马拉雅地区的三维速度模型，并利用高性能计算设施进行计算。研究首先使用2013年M_w5.8巴达瓦地震的记录数据对模型进行验证，然后针对1905年康格拉地震建立了多尺度震源模型，通过五层子源结构实现了0.001-20 Hz的宽频带模拟。

研究区域和地震特征方面，研究区域涵盖北纬28.0°至34.5°，东经73.5°至80.0°，包括四个构造带，从康格拉和昌巴向东延伸至贝阿斯河谷，构成康格拉-昌巴地震区。构造带包括昌巴推覆体、高喜马拉雅（HH）、小喜马拉雅的潘贾尔叠瓦带（PIZ）和亚喜马拉雅，从南到北依次排列。小喜马拉雅的PIZ和亚喜马拉雅的第三纪岩石被主边界逆冲断层（MBT）构造分离。

地震构造框架显示，印度-恒河平原的冲积层和亚喜马拉雅地区的褶皱第三纪被喜马拉雅褶皱前缘（HFF）构造分离。MCT由潘贾尔逆冲断层（PT）代表，它将推覆序列推到小喜马拉雅PIZ之上，并定义了昌巴推覆体的南部边界。小喜马拉雅构造的逆冲楔位于约10公里厚的主要极低级变沉积岩至沉积岩序列之上，构成昌巴推覆体。沿PT，推覆序列被迫覆盖PIZ，然后沿MBT逆冲到亚喜马拉雅的第三纪之上。PT和MBT之间的距离很小，只有几百米到几公里。在这种构造情况下，PT和MBT似乎共享向南的倾覆和构造运动。CN以东的喜马拉雅地区与1905年康格拉地震袭击的地区具有不同的地质背景。虽然CN和克什米尔序列（KS）在 stratigraphically 与特提斯喜马拉雅相当，位于HH带和小喜马拉雅PIZ的狭窄带之间，但在昌巴和克什米尔地区，HH结晶直接位于小喜马拉雅带之上并与之接触。MCT在地表表现为PT，而HH结晶和MCT在CN的东部闭合处 underlie CN。与Garhwal和Kumaun相比，康格拉地震地区的亚喜马拉雅要宽得多。这意味着在Garhwal-Kumaun地区沿MBT的亚喜马拉雅构造的底冲和在康格拉-昌巴地区沿PT的小喜马拉雅构造的底冲是汇聚的主要驱动因素。

介质建模方面，研究区域包括地形变化，从恒河平原的松散沉积物到喜马拉雅的坚硬岩石地形，因此必须纳入材料特性的横向变化。采用了深度达40公里的三维速度模型。扩展和改进了先前的三维速度模型，通过应用增强的空间插值和深度调整程序覆盖更广的研究区域。由于前模型的空问范围限于较小区域，我们通过执行空间插值将其覆盖范围扩展到我们的更大研究区域，并辅以参考的一维模型以确保扩展域中速度结构的一致性。实施了自适应插值策略——在平滑变化区使用线性插值，在地质复杂区使用三次插值——以确保连续的速度梯度，同时最小化伪影。深度调整依据更新的地球物理数据集和地层剖面，以保持真实的垂直对比，这对准确的地震波传播至关重要。此外，在印度-恒河平原内纳入了低速剪切波层（约475 m/s）至2公里深度，以更好地捕捉盆地放大效应。该地区的地形数据来自先进陆地观测卫星全球数字表面模型（AW3D30），提供30米的令人印象深刻的水平分辨率。该地区具有显著的高程变化，记录的最高高程为5152米。这种详细的地形信息对于准确模拟地震波放大至关重要，因为地震波与地形特征的相互作用可以深刻影响地震动特征。

研究方法上，使用谱元法（SEM）进行区域建模，该方法与使用形状函数的传统有限元法（FEM）有相似之处。在SEM中，采用高次拉格朗日插值定义位移场的基函数。具体地，每个元素由n次高斯-洛巴托-勒让德（GLL）控制点定义，包括端点。SEM的一个显著特征是质量矩阵的对角性。这种对角化源于使用GLL求积的插值点与分析一致。然而，需要注意的是，这种效率是以精度为代价的，与传统的FEM方法相比，精度可能略有妥协。SEM中的对角质量矩阵促进了高效的显式时间积分方案，使其特别适用于大规模动态模拟。尽管在精度上有所权衡，谱元法因其计算效率和用相对较少的元素处理复杂几何形状的能力而受到高度重视。

计算区域被视为三维弹性半空间，定义为满足超平面不等式一侧的所有点（即z≥0），假设材料域在平面之外无限延伸。在这里， traction-free 边界条件与吸收边界条件结合，成功模拟了无限或半无限空间的效果。半空间的表面，即地球表面，被认为是 traction-free，即σ·n=0，其中σ是应力张量，n是表面的单位法向量。计算域的边界可能导致波的人工反射，因此使用吸收边界条件来产生无限或半无限域的效果，特别是当要模拟波的传播时。这防止了波在边界处人工反射回计算域。这里使用的类型是完全匹配层边界条件，其中进入其中的波指数衰减并被阻止反射回系统。

在本研究中，开发的速度模型在SPECFEM Cartesian中实施，这是一组FORTRAN子程序，可用于地震动模拟。速度模型被离散化为六面体元素，这些元素使用雅可比矩阵映射到参考立方体上。如前所述，对于n阶多项式，需要n+1个高斯-洛巴托-勒让德（GLL）点。来自每个局部元素的贡献被组装形成全局矩阵，其中常微分方程控制全局系统的时间依赖性，表示为：Mü + C? + KU = F，其中U包含全局网格中所有网格点的位移向量，称为全局自由度。M是全局质量矩阵，C是全局吸收边界矩阵，K是全局刚度矩阵，F是源项。此外，方程的实践推进采用二阶有限差分格式。选择时间步长使得 courant number 小于0.3以确保系统稳定性。SPECFEM cartesian 包在AQUA集群中实施，该集群是印度理工学院马德拉斯高性能计算环境（HPCE）的一部分。该集群包括HPE Apollo XL170rGen10服务器，配备双英特尔至强金牌6248 20核2.5 GHz处理器，每个节点总共192 GB RAM。研究区域位于第43通用横轴墨卡托（UTM）带，覆盖面积为731.25km × 617.5km。该区域沿x和y方向离散为384个谱元素。沿表面每个方向使用8个消息传递接口（MPI），总共64个MPI处理器。沿深度进行元素细分，如前所述，有6个接口，进一步离散为层。按升序网格划分确保捕捉表面层的材料特性的微小变化，而底部的较粗网格有助于减少计算时间，因为特性变化不太显著。最小元素尺寸为0.107公里，采用的时间步长为0.002秒。网格生成和数据库所需时间分别为30分钟和12小时。

使用SPECFEM3D Cartesian代码（Komatitsch and Tromp, 1999）在计算域中执行数值模拟，计算域尺寸为731×618×40公里，分别沿x、y和z方向离散为384×384×80个六面体元素。表面的元素尺寸为1.9×1.6×0.1公里。为了优化计算效率，通过定义六个接口在垂直（z）方向应用"网格加倍"。这种方法通过在不牺牲显著模拟精度的情况下在更大深度使用更粗的网格来降低计算成本。虽然均匀网格将在深度提供更多网格点，但也会导致计算成本大幅增加，而模拟精度改善最小。因此，网格设计为在表面更细，那里材料特性的变化更明显。垂直域包括六个地层接口，从而定义了六个离散层。实施了分级网格方法，将最高空间分辨率应用于最上层，随着深度逐渐变粗。接口的位置主要对应于剪切波速度显著增加的深度。具体地，第6层（最顶层）以0.1公里间隔离散，第5层以0.2公里，第4层以0.4公里，三个基底层分别以1.0公里间隔离散。关于网格配置的更多技术细节在表1中提供。

表面高程使用来自先进陆地观测卫星全球数字表面模型（AW3D30）的地形数据纳入模型，该模型提供30米的水平分辨率。该地区表现出 substantial 高程变异性，最高点达到5152米。数字高程模型（DEM）在每个谱元素节点处插值，并分配自由边界条件。图3说明了研究区域的计算模型。模型中的总自由度为2.9亿。为确保数值稳定性，基于Courant稳定性准则选择0.002秒的时间步长。网格分辨率 tailored 以解析模型中的最低剪切波速度，使得能够模拟频率高达0.5 Hz的地震动时程。模拟使用印度理工学院马德拉斯的高性能计算设施进行。生成模型的合成地震动以分析研究区域内的地震情景。

在本研究中，对两次地震进行了地震动模拟：（1）2013年巴达瓦地震，在第4节中解释；（2）1905年康格拉地震，在第5节中解释。

为了验证三维计算模型（图3），生成了2013年5月1日12时27分（IST）发生的2013年巴达瓦地震（M_w5.8）的地震动，该地震袭击了巴达瓦克什米尔和喜马偕尔邦的边境地区。

震源模型方面，地震震中位于32.8°N, 76.3°E，位于巴达瓦东北23公里处。地震是浅地壳事件，震源深度23公里，持续约40秒。因此，本研究模拟了一个足够大的区域，以纳入大多数记录台站，这些台站与源-场地距离从74公里到332公里不等。事件的震源建模为有限断层，断层尺寸为10公里×7公里。断层机制确定为逆冲断层，走向角299°，滑动角80°，倾角24°，地震矩为3.6×10¹⁸ Nm，断层深度23公里。由于该地震的滑动分布尚未建立，采用了基于现实断层框架的随机滑动分布。这种空间随机场模型为生成旨在强地震动预测的情景地震提供了宝贵基础。使用 established 缩放关系估计震源尺寸，同时使用 spectral synthesis 方法生成空间可变的滑动分布。重要的是，每个异质滑动的实现都被约束以保持有限的静态应变能，确保应力保持有界。从这些滑动分布导出的强震动合成——假设简化的运动学破裂参数如恒定破裂速度和上升时间——突出了模型用于 robust 地震动预测的能力。

地震动方面，将所选台站的模拟速度时程与记录的地震动记录进行比较。记录加速度时程（对于场地类别B的查谟、场地类别C的穆克里安和场地类别C的鲁尔基）进行数值积分，获得速度时程，并进一步在0.05至0.8 Hz之间滤波，因为模拟地震动无法捕捉高频。图4显示时程同相并有效捕捉了每个事件所有方向的峰值。此外，使用傅里叶振幅谱（FAS）图比较了记录和模拟地震动的频率内容，如图4所示。记录和模拟数据的时间历史和频率振幅谱在很大程度上一致，两者显示出相似的整体趋势。尽管穆克里安的Joyner-Boore距离更大，但可以观察到查谟和穆克里安台站的速度时程振幅相当。这种现象主要归因于场地效应，特别是穆克里安相对于查谟显著更低的剪切波速度（Vs）。穆克里安降低的Vs导致由于较软、较不 rigid 的近地表沉积物中增加的场地共振和波 trapping 而增强的地震波放大。因此，局部场地放大机制抵消了更大距离衰减的影响，导致两个台站记录到相似的速度振幅。这表明模型在再现事件的时间演化和频率内容方面 reasonably 准确。因此，开发的计算模型可以进一步用于模拟1905年康格拉地震。

记录和模拟地震动数据之间的傅里叶振幅谱残差量化使用统计分析，如图5所示。估计是针对观测和模拟时程的对数缩放傅里叶振幅谱（FAS）获得的。每个频率的 misfit 或残差计算为记录和模拟振幅谱之间的差异。给定离散地面运动记录y(t)以Δt采样，计算其复傅里叶变换。对于每个分量，即EW、NS、垂直，计算残差谱。由于残差可以改变符号并跨越数量级，分析使用其绝对值的对数，这避免了对于y=0的未定义值。频率范围划分为0.5 Hz箱f_k，中心为f_k。在每个箱内，对数残差的中心趋势和 spread 为μ_k = (1/N_k) Σ_f∈Bk log|R(f)| 和 σ_k = sqrt((1/N_k) Σ_f∈Bk [log|R(f)| - μ_k]²)，其中N_k是落在B_k中的离散频率点的数量。绘制μ_k作为f_k处的散点，并添加长度为2σ_k的垂直段传达了平均偏差及其变异性。因此，这种紧凑表示暴露了数值模型在每个分量的频率带上 systematic 过度或不足预测模式。

1905年康格拉地震的确切震中由于缺乏地震学数据无法精确定位。然而，基于在康格拉-达兰萨拉地区观察到的破坏，文献中有 several 估计，震中的最佳估计约为32.18°N, 76.43°E。对事件的震级和震源深度进行了 several 推断，如'引言'中讨论。在本研究中，事件被认为是M_w7.8，震源深度18公里。基于各种研究，特别是Malik等人（2015），康格拉谷断层已被确定为该地区地震活动的主要来源。因此，本研究采用有限断层模型，断层尺寸为130公里×30公里。断层机制 characterized 为低角度逆冲断层，走向角299°，倾角19°，滑动角58°。

有限断层震源反演强调了跨断层平面的地震滑动的空间复杂性。表征这种复杂性的随机方法将滑动分布视为空间随机场，与已发布的有限源破裂模型一致。观测表明，von Karman自相关函数提供了 best 表示，相关长度与震源尺寸成比例，其中各向同性介质中使用极坐标的功率谱密度定义为P(k)∝(1+k²a²)^-H-1，其中H是Hurst指数。对于这项研究，历史地震的滑动模型尚未建立。因此，基于现实断层框架生成了随机滑动分布模型。开发了有限断层模型并划分为子源，关键震源参数——矩张量、破裂时间、上升时间和深度——为每个子源计算。最初，为指定的源函数和源几何创建了位错模型。然后使用 spectral synthesis 方法确定断层表面的滑动，其中滑动振幅谱通过空间自相关函数定义。每个子源的矩张量使用M₀=μ×D×A计算，能够为每个子源计算所有六个矩张量分量。对于此分析，土壤的剪切模量取为3.2×10¹⁰ N/m²。采用的滑动模型及其地理位置如图6所示，图7显示了半持续时间和时间偏移的分布。

这种空间随机场模型为生成旨在强地震动预测的情景地震提供了宝贵框架。可以使用 established 缩放关系估计震源尺寸，同时使用 spectral synthesis 方法生成空间可变的滑动分布。重要的是，每个异质滑动的实现都被约束以保持有限的静态应变能，确保应力保持有界。从这些滑动分布导出的强震动合成，假设简化的运动学破裂参数（例如，恒定破裂速度和上升时间），证明了模型用于 robust 地震动预测的潜力。速度时程（在0.001-0.5 Hz带通滤波器范围内滤波）为171个台站模拟，这些台站 within 所选区域有记录的宏观烈度测量值。

低频地震动模拟方面，模拟的速度时程按台站与断层的距离升序组织，提供了 seismic 波传播 across 不同 proximity 的结构化表示，以及 respective 傅里叶振幅谱（图8）。可以观察到获得的最大可解析频率为0.5 Hz。为清晰和重点，从总共171个台站中选择了12个，确保数据不仅具有代表性而且可理解。这种选择突出了 seismic 波形随着它们穿越与断层增加距离的空间演化。如图9所示，可以观察到P波到达时间的 discernible 延长随着距离增加，强调了 seismic 能量传输 over 扩展路径的 dispersive 性质。这种 dispersion 表现在波形的逐渐拉伸，这是 seismic 波与地壳异质结构相互作用的固有特征。并行地，seismic 波的振幅随着距离逐渐减小，这种现象归因于波传播过程中的能量衰减。这种衰减与理论预期一致，其中 seismic 波由于散射、吸收和几何扩展而损失能量。 together，这些观察提供了 seismic 波 across 各种地质介质行为的能量 dispersion 和衰减之间 intricate 平衡的清晰可视化。

宽频带地震动模拟的重新定义研究区域方面，前一节表明使用介质模型和震源模型模拟的时程仅对较小频率范围（<1 Hz）有效。然而，宽频带地震动对地震学和地震工程目的很重要。这突出了改进介质和震源模型以有效捕捉更高频率分量的必要性。由于高频地震动具有短波长，它需要介质和震源模型的精细离散化。精细震源模型有助于生成 realistic 高频内容，而精细介质模型确保这些频率传播而没有伪反射或数值 dispersion。为 address 此限制，提出了重新定义的研究区域用于模拟宽频带地震动，特别是近场台站。修订的域尺寸显著减小，现在涵盖北纬31.5°至33.0°，东经75.5°至77.5°，深度40公里。

数值模拟使用222×143×40公里的计算域进行，分别沿x、y和z方向离散为384×384×80个六面体元素。在表面，元素尺寸为0.57×0.37×0.1公里，确保关键区域的高空间分辨率。为优化计算效率，沿垂直（z）方向实施"网格加倍"策略。定义了六个接口，使得能够 progressively 更粗的网格随着深度，同时保持 near 表面更细的分辨率。这种方法显著降低计算成本而不妥协模拟精度，因为更深区域通常需要更少细节。网格的额外技术规格在表2中概述。

表面地形使用来自先进陆地观测卫星全球数字表面模型（AW3D30）的数据纳入模型。数字高程模型（DEM）在每个谱元素节点处插值，并应用自由边界条件以表示 realistic 表面相互作用。该地区的计算模型，包括高程细节，如图10所示。该模型包括2.9亿自由度，要求0.002秒的时间步长以确保数值稳定性，如由Courant稳定性准则所规定。网格设计为解析最低剪切波速度，使得能够模拟频率高达20 Hz的地震动时程。模拟使用印度理工学院马德拉斯的高性能计算设施进行，其中生成合成宽频带地震动以分析整个研究区域的地震情景。在开发重新定义研究的计算模型后，多尺度震源模型从基本有限断层模型导出。

多尺度震源模型方面，在本研究中，考虑矩形断层模型，具有矩形子源，尺寸Δl_i×Δw_i，如前一节讨论。前一节表明简单基本震源模型不足以捕捉高频地震动。为生成高频地震动，需要子源 reasonably 短的上升时间。在本研究中，多尺度震源模型基于假设开发，即几个 superimposed 层的相同破裂过程，每层均匀划分为子源，其尺寸和上升时间 progressively 减小可以生成高频 seismic 波。第一层中子源的上升时间与指定的上升时间分布对齐，保留定义源过程的关键特征。尺寸为R的子源的转角频率给出为f_c=β/(2πR)，其中β是剪切波速度。由于上升时间是转角频率的倒数，因此，尺寸为R的子源的上升时间给出为τ=2πR/β。从方程（7）明显看出，子源的上升时间与其尺寸成比例，因此不同层中子源的上升时间计算为：τ_i,j=C_i×τ_j，其中τ_i,j