In general, it will be a distributed RL network, with an overall \$|Z| \sim \sqrt{F}\$ characteristic, at least over the range of interest. This is due to the skin effect, a diffusion mechanism, which therefore has this impedance slope, or equivalently, a phase of 45°, or equal real and imaginary (reactive) components.
That fixes the slope of the function, but calculating the offset is nontrivial, and most generally will require a field simulation. Whether enough assumptions and approximations can be made to justify such a value, depends.
For practical purposes, I would suggest starting with an empirical or semi-empirical model. For example, laminated iron cores can be modeled with the generalized Steinmetz equation, \$p = C F^\alpha B^\beta\$ for some constants C, α and β, where p is the power density (usually mass density for transformer iron, volumetric (specific) density for high-frequency materials). These parameters are usually tabulated by the manufacturer and given in the datasheet.
Given core parameters, one can calculate the loss tangent of magnetizing current, or whatever equivalent core loss resistance, as reflected at the winding of interest, and extrapolate from there assuming skin effect is known and in play at the given frequency/range. (Preferably/hopefully, the loss equation is valid over the frequency range of interest, and so the loss resistance can be solved in terms of it, without having to make further assumptions at all. Mind that values may not be fit over a broad range, but only those relevant to mains frequencies, 50-400Hz for example.)
You may find the discussion here of interest,
Discussion | Calculators: Coilcraft SPICE Model Converter | Seven Transistor Labs
in many real cases (core loss and wire skin effect), we can approximate frequency-dependent resistance as a Warburg (diffusion) element, which in turn can be expressed with a circuit such as:
where the R and L values are in a geometric series.
Another example can be seen in this equivalent model,
Source: my site, CurveFit3.png
L5-L7 and R4, R6, R7 approximate the shallow-rising impedance in the 20-400 kHz range. The part uses a very thin strip-type core material (nanocrystalline), which therefore exhibits skin effect, which dominates its impedance in this range.
As for flux:
Skin effect is a self-shielding effect, preventing magnetic field from reaching the center of the conductor. Cross-sectional area therefore drops at a corresponding rate (hence this curve manifests for both conductor and core losses). This might cause less flux to be seen by the core overall, or the remaining cross-section is driven more strongly (i.e., the same given magnetic field strength, but \$A_e\$ is smaller), and thus driven closer to saturation, which acts to further increase core losses.
For solid materials, I'm not aware of a simple calculation, even for common geometries such as a cylinder inside a solenoid coil. Further reading might be found in texts on induction heating; though I don't have any recommendations handy unfortunately.
