A random parameter field generator within FEFLOW is planned for the future. It has not been developed yet.
Accordingly, you need to use an external generator like the one you indicated. You could use the output of the generator as input in FEFLOW to parametrize your model.
One option is to convert the generator output into the ASCII based dat-file format readable by FEFLOW. This file can be used as a map for parameter linking (regionalization). The disadvantage is the manual work. I would expect that this task is tedious and time-consuming particularly if you have to carry out a large number of simulations.
Another possibility is the development of a script (Python) or plugin (C/C++). I suggest to use Python, where Kernel Control Methods are provided. Generally, Kernel control methods enable the user to start/stop the kernel, to load/save documents, and to control the simulation. You could use kernel control method [font=courier][color=blue]doc.startSimulator()[/color][/font] to start the simulator within a loop ranging from 0 to n_number_of_simulations-1. Before you start the simulator you could assign the random parameter field to the mesh elements
In this context, please note that available regionalization methods are not available for Python. Instead you need to assign the K-values via
[font=courier][color=blue]
extern void IfmSetMatXConductivityValue3D (IfmHandle, int, double);
extern void IfmSetMatYConductivityValue3D (IfmHandle, int, double);
extern void IfmSetMatZConductivityValue3D (IfmHandle, int, double);
[/color][/font]
You could derive the K-values of each element by taking advantage of using the FEFLOW API and GIS API (e.g. arcpy) within a single script. Since the mesh do not change over time (I assume you do not use the free & moveable method) you need to export the element centroids as a point shapefile from the Data panel only one time (not from the script).
Within each loop you may convert the parameter field into a GIS raster file or TIN (Triangulated Irregular Network). To transfer (project) the K-values of the raster / TIN to the shapefile attributes you could use the arcpy function [font=courier][color=blue]AddSurfaceInformation_3d (in_feature_class, in_surface, out_property, {method}, {sample_distance}, {z_factor}, {pyramid_level_resolution}, {noise_filtering})[/color][/font]. In former times, the function was also called Surface Spot.
After the K-values are written to the input feature's attribute table you may use the FEFLOW API function [font=courier][color=blue]IfmFindElementAtXYZ (IfmHandle, double x, double y, double z);[/color][/font] to derive the element ID's which correspond to the points. As soon you derived the element ID you may finally assign the K-value based on the three assignment functions I have shown above.
In general the workflow seems a bit complicated, but the script would only have a few lines of code.